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This article concerns numerical simulations of the dynamics of particles immersed in a 
continuum solvent. As prototypical systems, we consider colloidal dispersions of spherical 
particles and solutions of uncharged polymers. After a brief explanation of the concept of 
hydrodynamic interactions, we give a general overview of the various simulation methods that 
have been developed to cope with the resulting computational problems. We then focus on the 
approach we have developed, which couples a system of particles to a lattice-Boltzmann model 
representing the solvent degrees of freedom. The standard D3Q19 lattice-Boltzmann model is 
derived and explained in depth, followed by a detailed discussion of complementary methods 
for the coupling of solvent and solute. Colloidal dispersions are best described in terms of 
extended particles with appropriate boundary conditions at the surfaces, while particles with 
internal degrees of freedom are easier to simulate as an arrangement of mass points with 
frictional coupling to the solvent. In both cases, particular care has been taken to simulate 
thermal fluctuations in a consistent way. The usefulness of this methodology is illustrated by 
studies from our own research, where the dynamics of colloidal and polymeric systems has 
been investigated in both equilibrium and nonequilibrium situations. 
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1 

Introduction 

The term "soft condensed matter" generally refers to materials which possess addi- 
tional "mesoscopic" length scales between the atomic and the macroscopic scales [1- 
8]. While simple fluids are characterized by the atomic size (3 x 10~ 10 m), soft- 
matter systems contain one or more additional length scales, typically of order 
10~ 9 m — 10~ 6 TO. There are many examples of matter with mesoscale structure, 
including suspensions, gels, foams, and emulsions; all of these are characterized by 
viscoelastic behavior, which means a response that is fluid-like on long time scales 
but solid-like on shorter time scales. The two prototype systems considered in this 
article are colloidal dispersions of hard particles, where the additional length scale is 
provided by the particle size, and polymer systems, where the length scale is the size 
of the macromolecule. The main difference between these two systems is the pres- 
ence of internal degrees of freedom in the polymer, such that a statistical description 
is necessary even on the single-molecule level. Many additional soft-matter sys- 
tems exist. For example, dispersions may not only contain spherical particles, but 
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rather rod-like or disk-like objects. For polymers, there are many possible molec- 
ular architectures; in addition to simple linear chains, there are rings, stars, combs, 
bottle-brush polymers and dendrimers. Polymers may also self-assemble into two- 
dimensional membranes, either free or tethered, which are of paramount biological 
importance in cell membranes, vesicles, and red blood cells. 

Strongly non-linear rheology is characteristic of soft matter. In simple fluids, it is 
difficult to observe any deviations from Newtonian behavior, which is well described 
by the hydrodynamic equations of motion with linear transport coefficients that de- 
pend only on the thermodynamic state. Indeed, Molecular Dynamics simulations [9] 
have revealed that a hydrodynamic description is valid down to astonishingly small 
scales, of the order of a few collisions of an individual molecule. This means that one 
would have to probe the system with very short wave lengths and very high frequen- 
cies, which are typically not accessible to standard experiments (with the exception 
of neutron scattering [10]), and even less in everyday life. However, in soft-matter 
systems microstructural components (particles and polymers for example) induce 
responses that depend very much on frequency and length scale. These systems are 
often referred to as "complex fluids". 

The nonlinear rheological properties of soft matter pose a substantial challenge 
for theory [11]. Therefore, the study of simple model systems is often the only way 
to make systematic progress. Numerical simulations allow us to follow the dynamics 
of model systems without invoking the uncontrolled approximations that are usu- 
ally required by purely analytical methods. Simulations can be used to isolate and 
investigate the influence of microstructure, composition, external perturbation and 
geometry in ways that cannot always be duplicated in the laboratory. In particular, 
simulations can provide a well-defined test bed for theoretical ideas, allowing them 
to be evaluated in a simpler and more rigorous environment than is possible ex- 
perimentally. Finally they can provide more detailed and direct information on the 
particle dynamics and structure than is typically possible with experimental mea- 
surements. 

In this article we will focus on systems which comprise particles, with or without 
internal degrees of freedom, suspended in a simple fluid. We will first outline the nec- 
essary ingredients for a theoretical description of the dynamics, and in particular ex- 
plain the concept of hydrodynamic interactions (HI). Starting from this background, 
we will provide a brief overview of the various simulation approaches that have been 
developed to treat such systems. All of these methods are based upon a description 
of the solute in terms of particles, while the solvent is taken into account by a simple 
(but sufficient) model, making use of the fact that it can be described as a Newtonian 
fluid. Such methods are often referred to as "mesoscopic". We will then describe 
and derive in some detail the algorithms that have been developed by us to couple a 
particulate system to a lattice-Boltzmann fluid. The usefulness of these methods will 
then be demonstrated by applications to colloidal dispersions and polymer solutions. 
Some of the material presented here is a summary of previously published work. 
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Particle-fluid systems 



2.1 



Coarse-grained models 

The first step towards understanding systems of particles suspended in a solvent is 
the notion of scale separation. Colloidal particles are much larger (up to lCT 6 ™) 
than solvent molecules, and their relaxation time (up to 10s), given by the time the 
particle needs to diffuse its own size, is several orders of magnitude larger than the 
corresponding solvent time scale (10~ 12 s) as well. A similar separation of scales 
holds for solute with internal degrees of freedom, like polymer chains, membranes, 
or vesicles, as far as their diffusive motion or their global conformational reorgani- 
zation is concerned. However, these systems contain a hierarchy of length and time 
scales, which can be viewed as a spectrum of internal modes with a given wavelength 
and relaxation time. For the long-wavelength part of this spectrum, the same scale 
separation holds, but when the wavelength is comparable with the solvent size molec- 
ular interactions become important. However, on scales of interest, these details can 
usually be lumped into a few parameters. The solute is then modeled as a system of 
"beads" interacting with each other via an effective potential. The beads should be 
viewed as collections of atomic-scale constituents, either individual atoms or func- 
tional groups, which have been combined into a single effective unit in a process 
known as "coarse-graining". On the scale of the beads, the solvent may be viewed 
as a hydrodynamic continuum, characterized by its shear viscosity and temperature. 
The flow of soft matter is usually isothermal, incompressible, and inertia-free (zero 
Reynolds number). Then, the most natural parameter to describe solute-solvent cou- 
pling in polymeric systems is the Stokes friction coefficient of the individual beads, 
or (in the case of anisotropic subunits) the corresponding tensorial generalization. 
The effective friction coefficient is the lumped result of a more detailed, or "fine- 
grained", description at the molecular scale, as is the bead-bead potential, which 
should be viewed as a potential of mean force. In order to make contact with experi- 
mental systems, these parameters must be calculated from more microscopic theories 
or simulations, or deduced from experimental data. 

Further simplification may arise from the type of scientific question being ad- 
dressed by the coarse-grained model. If the interest is not in specific material prop- 
erties of a given chemical species, but rather in generic behavior and mechanisms, as 
is the case for the examples we will discuss in this article, then the details of the pa- 
rameterization are less important than a model that is both conceptually simple and 
computationally efficient. The standard model for particulate suspensions is a system 
of hard spheres, while for polymer chains the Kremer-Grest model [12] has proved 
to be a valuable and versatile tool. Here, the beads interact via a purely repulsive 
Lennard-Jones (or WCA [13]) potential, 
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while consecutive beads along the chain are connected via a FENE ("finitely exten- 
sible nonlinear elastic") potential, 



with typical parameters Rq/u = 1.5, ka 2 /e = 30. Chain stiffness can be incor- 
porated by an additional bond-bending potential, while more complicated architec- 
tures (like stars and tethered membranes) require additional connectivity. Poor sol- 
vent quality is modeled by adding an attractive part to the non-bonded interaction, 
while Coulomb interactions (however without local solvent polarization) can be in- 
cluded by using charged beads. 

There are two further ingredients which any good soft-matter model should take 
into account: On the one hand, thermal fluctuations are needed in order to drive 
Brownian motion and internal reorganization of the conformational degrees of free- 
dom. There are however special situations where thermal noise can be disregarded, 
and ideally the simulation method should be flexible enough to be able to turn the 
noise both on and off. On the other hand, hydrodynamic interactions, which will be 
the subject of the next subsection, need to be taken into account in most circum- 
stances. 



Hydrodynamic interactions 

The term "hydrodynamic interactions" describes the dynamic correlations between 
the particles, induced by diffusive momentum transport through the solvent. The 
physical picture is the same, whether the particle motion is Brownian {i.e. driven by 
thermal noise) or the result of an external force {e.g. sedimentation or electrophore- 
sis). The motion of particle i perturbs the surrounding solvent, and generates a flow. 
This signal spreads out diffusively, at a rate governed by the kinematic viscosity of 
the fluid rjkin = v/P (V is the solvent shear viscosity and p is its mass density). On 
interesting (long) time scales, only the transverse hydrodynamic modes [14] remain, 
and the fluid may be considered as incompressible. The viscous momentum field 
around a particle diffuses much faster than the particle itself, so that the Schmidt 
number 



is large. In a molecular fluid, D is the diffusion coefficient of the solvent molecules 
and Sc <~ 10 2 — 10 3 , while for soft matter, where D is the diffusion coefficient of the 
polymer or colloid, Sc is much larger, up to 10 6 for micron size colloids. The con- 
dition Sc ^> 1 is an important restriction on the dynamics, which good mesoscopic 
models should satisfy. Consider two beads i and j, separated by a distance r^. The 
momentum generated by the motion of i with respect to the surrounding fluid reaches 
bead j after a "retardation time" t ~ rfj/r/kin- After this time the motion of j be- 
comes correlated with that of bead i. However, during the retardation time the beads 
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have traveled a small distance <~ V Dt ~ r^/VSc, which is negligible in compari- 
son with rij . Therefore, it is quite reasonable to describe the Brownian motion of the 
beads neglecting retardation effects, and consider their random displacements to be 
instantaneously correlated. 

These general considerations suggest a Langevin description (stochastic differ- 
ential equation) for the the time evolution of the bead positions rf. 

J t T ia = Vic,j0 F jP + A ia , (4) 
3 

where a, (3 are Cartesian indexes and the Einstein summation convention has been 
assumed. F^ is the conservative force acting on the jth bead, while the mobility 
tensor describes the velocity response of particle i to FJ; both F^ and ^ depend 
on the configuration of the N particles, r N . The random displacements (per unit 
time) A ia are Gaussian white noise variables [15, 16] satisfying the fluctuation- 
dissipation relation 

(Aa) = 0, (5) 

(A ia {t)A jf3 (t')) = 2k B T^ a .^S(t - t'); (6) 

here T is the temperature and k B denotes the Boltzmann constant. In general, the mo- 
bility tensor, ^iajp depends on particle position, in which case the stochastic integral 
should be evaluated according to the Stratonovich calculus. This is because during a 
finite time step a particle samples different mobilities as its position changes. Numer- 
ically simpler is the Ito calculus [15, 16], which uses the mobility at the beginning 
of the time step. In this case an additional term needs to be added to Eq. 4, 

|r fa = ]T Wa^Jft + k B Tj2 + A a . (7) 

3 3 • 7/3 

In cases where the divergence of the mobility vanishes, d^iajp/drjp = 0, the Ito 
and Stratonovich interpretations coincide. 

Processes generated by Eq. 4 (Stratonovich) or Eq. 7 (Ito) sample trajectories 
from a probability distribution in the configuration space of the beads, P (r N , t) , 
which evolves according to a Fokker-Planck equation (Kirkwood diffusion equa- 
tion): 

^P(r N ,t)=£P(r N ,t), (8) 
with the Fokker-Planck operator 

£ = E^«(^T^--^). (9) 

Writing the forces as the derivative of a potential, 

Fr a = ~^V(r N ), (10) 

u 1 ia 
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it follows that the Boltzmann distribution, P oc cxp (—V/k B T), is the stationary 
solution of the Kirkwood diffusion equation; in other words the model satisfies the 
fluctuation-dissipation relation. 

The mobility tensor can be derived from Stokes-flow hydrodynamics. Consider 
a set of spherical particles, located at positions r ; with radius a, surrounded by a 
fluid with shear viscosity 77. Each of the particles has a velocity v», which, as a 
result of stick boundary conditions, is identical to the local fluid velocity on the 
particle surface. The resulting fluid motions generate hydrodynamic drag forces Ff , 
which at steady state are balanced by the conservative forces, + = 0. The 
commonly used approximation scheme is a systematic multipole expansion, similar 
to the analogous expansion in electrostatics [17-21]. For details, we refer the reader 
to the original literature [17], where the contributions from rotational motion of the 
beads are also considered. As a result of the linearity of Stokes flow, the particle 
velocities and drag forces are linearly related, 

Via = _ ViajpFjp = X! Via,j0 F 3f3- (H) 
3 3 

Since /i^ describes the velocity response of particle i to the force acting on particle 
j, it must be identical to the mobility tensor appearing in the Langevin equation. 

In general the mobility matrix is a function of all the particle coordinates, but to 
leading order, it is pairwise additive: 

djjl (1-%) ( vjjYjA (!-<%) a 2 / Jt^jA nr . 

where = — rj, and 1 denotes the unit tensor. The hydrodynamic interaction 
is long-ranged and therefore has a strong influence on the collective dynamics of 
suspensions and polymer solutions. This approximate form for the mobility matrix 
follows from the assumption that the force density on the sphere surface is constant. 
A point multipole expansion, by contrast, generates the Oseen (l/ry) interaction at 
lowest order [22], and can lead to non-positive-definite mobility matrices [23]. Thus 
the simplest practical form for the hydrodynamic interaction is the Rotne-Prager ten- 
sor [23] given in Eq. 12. Both the Oseen and Rotne-Prager mobilities are divergence 
free, and therefore there is no distinction between Ito and Stratonovich interpreta- 
tions. However, at higher orders in the multipole expansion, the divergence is non- 
zero [24]. 

2.3 

Computer simulation methods and models 

In this section we briefly summarize the Brownian dynamics algorithm and its close 
cousin Stokesian Dynamics. We then outline the motivation and development of sev- 
eral mesoscale methods, some of which are reviewed elsewhere in this series. 
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2.3.1 

Brownian Dynamics 

Brownian dynamics is conceptually the most straightforward approach [25, 26]. 
Starting from the Langevin equation for the particle coordinates, Eq. 4, and dis- 
cretizing the time into finite length steps h, gives a first-order (Euler) update for the 
particle positions, 

r ia (t + h) = r ia (t) + ^2 ViajpQjp- (13) 



Note that in Eq. 13 we have assumed a divergence-free mobility tensor 

Although the number of degrees of freedom has been minimized, this approach 
is computationally intensive, and imposes severe limitations on the size of the system 
that can be studied. Since every particle interacts with every other particle, the calcu- 
lation of the mobility matrix scales as 0(N 2 ), where N is the number of Brownian 
particles. In addition, the covariance matrix for the random displacements requires 
a Cholesky decomposition of the mobility matrix, which scales as 0(N 3 ) [27]. The 
computational costs of Brownian dynamics are so large that even today one cannot 
treat more than a few hundred Brownian particles [28]. 

"Stokesian Dynamics" [29] is an improved version of Brownian dynamics, in 
which the mobility tensor takes into account short-range (lubrication) contributions 
to the hydrodynamic forces. It also improves the far- field interactions by including 
contributions from torques and stresslets, although still higher moments are needed 
for accurate results in concentrated suspensions [19]. Stokesian Dynamics is even 
more computationally intensive than Brownian dynamics; the determination of the 
mobility tensor is already an 0(N 3 ) process. 

However, there have been two important improvements in efficiency. First, Fix- 
man [25, 26, 30, 31] has proposed an approximation to aia,jp by a truncated expan- 
sion in Chebyshev polynomials, which has a more favorable scaling than Cholesky 
decomposition. Second, the long-range hydrodynamic interactions can be calculated 
by Fast Fourier Transforms [32-37], or hierarchical multipole expansions [20]. Ac- 
celerated Brownian Dynamics and Stokesian Dynamics algorithms scale close to lin- 
early in the number of particles, and their full potential is not yet explored. However, 
it should be noted that all these methods are based upon an efficient evaluation of the 
Green's function for the Stokes flow, which depends on the global boundary condi- 
tions. For planar boundaries, solutions are available [21, 38-40], but a more general 
shape requires a numerical calculation of the Green's function between a tabulated 
set of source and receiver positions [41]. 



j 



Here qi a are random variables with 

(Qia) = 
{qi a qjp) = 8ij5 a p, 
while the matrix (Ti a ,ii3 satisfies the relation 



(14) 
(15) 




(16) 



k 
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2.3.2 

Mesoscale methods 

In view of the computational difficulties associated with Brownian Dynamics, several 
"mesoscale" methods have been developed recently. The central idea is to keep the 
solvent degrees of freedom, but to describe them in a simplified fashion, such that 
only the most salient features survive. As we have already seen, it is in principle 
sufficient to describe the solvent as a Navier-Stokes continuum, or by some suitable 
model which behaves like a Navier-Stokes continuum on sufficiently large length 
and time scales. At least asymptotically, the solvent dynamics must be described by 
the equations 

d t p + d a (pu a ) = 0, 

d t {pu a ) + dp {pu a u p ) + d a p = dp<j af3 + dpa^ + f a , (17) 

where p is the mass density, pu the momentum density, p the thermodynamic pres- 
sure, f an external force density applied to the fluid, a the viscous stress tensor, and 
a* the fluctuating (Langevin) stress [42], whose statistical properties will be dis- 
cussed in later chapters of this article. The viscous stresses are characterized by the 
shear and bulk viscosities, r\ and r] v , which we will assume to be constants, indepen- 
dent of thermodynamic state and flow conditions: 



The advantage of such approaches is their spatial locality, resulting in favorable 
O(N) scaling, combined with ease of implementation and parallelization. The disad- 
vantage is the introduction of additional degrees of freedom, and of additional (short) 
time scales which are not of direct interest. The coupling between solvent and solute 
varies from method to method. However, in all cases one takes the masses and the 
momenta of the solute particles explicitly into account, and makes sure that the total 
momentum is conserved. 

Lattice models (Navier-Stokes, lattice Boltzmann) simulate a discretized field 
theory in which thermal fluctuations can be added, but also avoided if desired. Par- 
ticle methods (Molecular Dynamics, Dissipative Particle Dynamics, Multi-Particle 
Collision Dynamics) simulate a system of interacting mass points, and therefore ther- 
mal fluctuations are always present. The particles may have size and structure or 
they may be just point particles. In the former case, the finite solvent size results 
in an additional potential of mean force between the beads. The solvent structure 
extends over unphysically large length scales, because the proper separation of scale 
between solute and solvent is not computationally realizable. In dynamic simulations 
of systems in thermal equilibrium [43], solvent structure requires that the system be 
equilibrated with the solvent in place, whereas for a structureless solvent the solute 
system can be equilibrated by itself, with substantial computational savings [43]. Fi- 
nally, lattice models have a (rigorously) known solvent viscosity, whereas for particle 
methods the existing analytical expressions are only approximations (which however 
usually work quite well). 
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These considerations suggest that lattice methods are somewhat more flexible 
and versatile for soft-matter simulations. On the other hand, the coupling between 
solvent and immersed particles is less straightforward than for a pure particle sys- 
tem. The coupling between solid particles and a lattice-based fluid model will be 
discussed in detail in Sec. 4. 

2.3.3 

Molecular Dynamics (MD) 

Molecular Dynamics is the most fundamental approach to soft-matter simulations. 
Here the solute particles are immersed in a bath of solvent molecules and Newton's 
equations of motion are solved numerically. In this case, it is impossible to make the 
solvent structureless - a structureless solvent would be an ideal gas of point parti- 
cles, which never reaches thermal equilibrium. Furthermore, the model interaction 
potentials are stiff and considerable simulation time is spent following the motion 
of the solvent particles in their local "cages". These disadvantages are so severe that 
nowadays MD is rarely applied to soft-matter systems of the type we are discussing 
in this article. 

2.3.4 

Dissipative Particle Dynamics (DPD) 

Dissipative Particle Dynamics, which has become quite popular in the soft-matter 
community [44-56], was developed to address the computational limitations of 
MD. A very soft interparticle potential, representing coarse-grained aggregates 
of molecules, enables a large time step to be used. Furthermore, a momentum- 
conserving Galilean-invariant thermostat is included, representing the degrees of 
freedom that have been lost in the coarse-graining process. Practically, these two 
parts are unrelated, such that it is legitimate to apply the DPD thermostat to a stan- 
dard MD system. The DPD thermostat is consistent with macroscopic isothermal 
thermodynamics. Since this already introduces interparticle collisions, it is possible 
to run DPD using an ideal gas solvent and still achieve thermal equilibrium.. 

The key innovation in DPD is to apply the thermostat to particle pairs. A fric- 
tional damping is applied to the relative velocities between each neighboring pair, 
and a corresponding random force is added in a pairwise fashion also, such that 
Newton's third law holds exactly. The implementation is as follows. We define two 
functions £(r) > 0, the relative friction coefficient for particle pairs with interpar- 
ticle distance r, and a(r) > 0, characterizing the strength of the stochastic force 
applied to the same particle pair. The fluctuation-dissipation theorem requires that 

a 2 (r) = k B TC(r). (19) 

The functions have compact support, so that only near neighbors need be taken into 
account. 

The frictional force on particle i is determined by projecting the relative veloci- 
ties onto the interparticle separation (fjj = r^/ jr^ |): 
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^ = -T,S(r ij )[(v i -v j )-r ij ]r ij , (20) 

3 

which conserves momentum exactly, J2i — 0- Similarly, the stochastic forces 
are directed along the interparticle separation, again so that momentum is conserved 
pair-by-pair, 

pf = I>( r y )%•(*)*«■ (2D 

3 

The noise r]ij satisfies the relations r]ij = r/ji, (rjij) = 0, and 

(Vij(t)Vki(t')) = 2(5ikSji + 5 u 5 jk )S(t - t'), (22) 

such that different pairs are statistically independent and ^\ F{ = 0. The equations 
of motion for a particle of mass rrn and momentum p i are: 

d 1 

~77 v i = Pi, (23) 

at m,i 

^p 4 = F 4 c + Ff + F{. (24) 

Exploiting the relation between this stochastic differential equation and its Fokker- 
Planck equation, it can be shown that the fluctuation-dissipation theorem holds [46], 
and that the method therefore simulates a canonical ensemble. DPD can be extended 
to thermalize the perpendicular component of the interparticle velocity as well, 
thereby allowing more control over the transport properties of the model [49, 57]. 



2.3.5 

Multi-Particle Collision Dynamics (MPCD) 

This method [58-62] works with a system of ideal-gas particles and therefore has 
no artificial depletion forces. Free streaming of the particles, 

r t (t + h) = Ti (t) +hvi(t), (25) 

alternates with momentum and energy conserving collisions, which are implemented 
via a Monte Carlo procedure: 

Sub-divide the simulation volume into a regular array of cells. 

• For each cell, determine the set of particles residing in it. For one particular cell, 
let these particles be numbered i = 1, . . . , n. Then in each box: 

• Determine the local center-of-mass velocity: 



1 n 

71 ^ ' 



vcm = - > Vj. (26) 



n 



• For each particle in the cell, perform a Galilean transformation into the local 
center-of-mass system: 

Vj = Vj - v C M- (27) 



12 



Burkhard Diinweg and Anthony J. C. Ladd 



Within the local center-of-mass system, rotate all velocities within the cell by a 
random rotation matrix R: 



By suitable random shifts of the cells relative to the fluid, it is possible to recover 
strict Galilean invariance [59, 60]. MPCD results in hydrodynamic behavior on 
large length and time scales, and is probably the simplest and most efficient parti- 
cle method to achieve this. 

2.3.6 

Lattice Boltzmann (LB) 

Here one solves the Boltzmann equation, known from the kinetic theory of gases, in 
a fully discretized fashion. Space is discretized into a regular array of lattice sites, 
time is discretized, and velocities are chosen such that one time step will connect 
only nearby lattice sites. Free streaming along the lattice links alternates with local 
on-site collisions. Care must be taken to restore isotropy and Galilean invariance 
in the hydrodynamic limit, and asymptotic analysis is an indispensable tool in this 
process. Further details will be provided in the following chapters. 

2.3.7 

Navier-Stokes 

It is possible to start from a discrete representation of Eqs. 17 but this has not been 
particularly popular in soft-matter simulations, due to the difficulty of including ther- 
mal fluctuations (but see Ref. [63]). Finite-difference methods share many techni- 
cal similarities with LB and are roughly comparable in terms of computational re- 
sources. However, to our knowledge, no detailed benchmark comparisons are avail- 
able as yet. In order to be competitive with LB, we believe that the solver must 
(i) make sure that mass and momentum are conserved within machine accuracy, as 
is the case for LB, and (ii) not work in the incompressible limit, in order to avoid 
the costly non-local constraints imposed by the typical Poisson solver for the pres- 
sure. The incompressible limit is an approximation, which eliminates the short time 
scales associated with wave-like motion. However, in soft matter the solute particles 
must be simulated on short inertial time scales, which requires that the solvent is 
simulated on rather short time scales as well. For this reason, we believe that enforc- 
ing an incompressibility constraint does not pose a real advantage, and it is instead 
preferable to allow for finite compressibility, such that one obtains an explicit and 
local algorithm. This idea is analogous to the Car-Parrinello method [64], where the 
Born-Oppenheimer constraint is also discarded, in favor of an approximate but ade- 
quate separation of time scales. For simulations of soft-matter systems coupled to a 
Navier-Stokes background, see Refs. [65-74]. 



= Rv*. 

Transform back into the laboratory system: 
v » = v- + v CM . 



(28) 



(29) 
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The fluctuating lattice-Boltzmann equation 

The motivation for the development of lattice gas cellular automata (LGCA) [75, 76] 
was to apply a highly simplified Molecular Dynamics (MD) to simulations of hydro- 
dynamic flows. In LGCA, particles move along the links of a regular lattice, typically 
cubic or triangular. Each lattice direction is encoded with a label i and a vector /iCj 
connects neighboring pairs of sites. During each time step h, all particles with a di- 
rection i are displaced /iCj to an adjacent lattice site; thus is the (constant) velocity 
of particles of type i. Interparticle interactions are reduced to collisions between par- 
ticles on the same lattice site, such that the conservation laws for mass and momen- 
tum are satisfied; in single speed LGCA models [75, 76], mass conservation implies 
energy conservation as well. The lattice-Boltzmann (LB) method [77-80] was de- 
veloped to reduce the thermal noise in LGCA, which requires extensive averaging to 
obtain statistically significant results. 

The LB model preserves the structural simplicity of LGCA, but substitutes an 
ensemble-averaged collision operator for the detailed microscopic dynamics of the 
LGCA. The hydrodynamic flow fields develop without thermal noise, but the under- 
lying connection with statistical mechanics is lost (Sec. 3.1). The LB model turns out 
to be more flexible than LGCA, and there is now a rich literature that includes ther- 
mal [8 1-86] and multiphase flows, involving both liquid-gas coexistence and mul- 
ticomponent mixtures [87-96]. In the present article, we will consider only single- 
phase flows of a single solvent species, such that we can describe the dynamics in 
terms of a single particle type. The algorithm can be summarized by the equation 

i>i(r + Cih,t + h) = v*{r,t) = v t {v,t) + A t (v(r,t)) , (30) 

where z^(r, t) is the number of particles that, at the discrete time t just prior to col- 
lision, reside at the lattice site r, and have velocity Cj; v*(r, t) indicates the velocity 
distribution immediately after collision. The difference A{ between the pre- and post- 
collision states is called the "collision operator" and depends on the complete set of 
populations at the site v(v.t). The left hand side of Eq. 30 describes the advection 
of the populations along the links connecting neighboring lattice sites. The velocity 
set Cj is chosen such that each new position r + Cj/i is again at a lattice site; Ci = 
is possible. 

For simplicity and computational efficiency, the number of velocities should be 
small. Therefore the set of velocities, c,, is typically limited to two or three neighbor 
shells, chosen to be compatible with the symmetry of the lattice. In two dimensions 
a single shell of 6 neighbors is sufficient for hydrodynamic flows, but a single set of 
cubic lattice vectors leads to anisotropic momentum diffusion, even at large spatial 
scales. Thus LB models employ a judicious mixture of neighboring shells, suitably 
weighted so that isotropy is recovered. We use the classification scheme introduced 
by Qian et al. [97]: for instance D2Q9 refers to an LB model on a square lattice 
in two dimensions, using nine velocities (zero, four nearest neighbors, four next- 
nearest neighbors), while D3Q19 indicates a three-dimensional model on a sim- 
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pie cubic lattice with nineteen velocities (zero, six nearest neighbors, twelve next- 
nearest neighbors). 

3.1 

Fluctuations 

The difference between lattice gas and lattice Boltzmann lies in the nature of the v±. 
In a lattice gas is a Boolean variable (i.e. only the values zero and one are allowed), 
while in the lattice-Boltzmann equation it is a positive real-valued variable. In Sec. 
3.6 we will consider the case where Vi is a large positive integer, a conceptual model 
we call a "Generalized Lattice Gas" (GLG). Thinking of these models as a simplified 
Molecular Dynamics, and considering fluctuations in v^, it becomes clear what the 
key difference between LGCA and the LB equation is. We define a dimensionless 
"Boltzmann number", Bo, by the fluctuations in v j at a single site, 



where (. . .) denotes the ensemble average. One could define Boltzmann numbers for 
other observables, but they would all produce similar values. The important point 
is that Bo tells us how coarse-grained the model is, compared to microscopic MD: 
Bo <~ 1 (the maximum value) corresponds to a fully microscopic model where fluc- 
tuations are of the same order as the mean. This is exactly the case for lattice gas 
cellular automata, which should therefore be viewed as a simplified, but not coarse- 
grained MD. Conversely, deterministic LB algorithms, at sufficiently small Reynolds 
numbers, and with time-independent driving forces, bring the system to a stationary 
state with well-defined values for the j/,. In other words, they are characterized by 
Bo = 0, which is the minimal value, corresponding to entirely deterministic physics. 

Originally, LGCA and LB algorithms were developed to simulate macroscopic 
hydrodynamics. Here, a large Boltzmann number (order 1) is undesirable, since the 
hydrodynamic behavior is only revealed after extensive sampling. For many macro- 
scopic applications a deterministic LB simulation at Bo = is hence entirely appro- 
priate. In reality, however, the Boltzmann number is finite, since the spatial domain 
in the physical system corresponding to a single lattice site is also finite. In soft- 
matter applications the spatial scales are so small that these fluctuations do need to 
be taken into account, although in many cases Bo is fairly small. This suggests it 
would be advantageous to introduce small thermal fluctuations into the LB algo- 
rithm, in a controlled fashion, by means of a stochastic collision operator [98-100]. 
The fluctuation-dissipation relation can be satisfied by enforcing consistency with 
fluctuating hydrodynamics [42] on large length and time scales. An important refine- 
ment is to thermalize the additional degrees of freedom that are not directly related to 
hydrodynamics [101], which leads to equipartition of fluctuation energy on all length 
scales. A comprehensive understanding of these approaches in terms of the statistical 
mechanics of LB systems has been achieved only recently [102]. 



Bo 




(31) 
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The number variables, v i7 can be connected to the hydrodynamic fields, mass 
density p(r, t), momentum density j(r, t), and fluid velocity u(r, t) (j = pu), by 
introducing the mass of an LB particle, m p , and the mass density parameter 

M = (32) 

here b is the lattice spacing and a three-dimensional lattice has been assumed. We 
then use the mass densities of the individual populations, 

rii(r,t) = fJ,Vi(r,t), (33) 

to re-write the LB equation as 

n t {r + Ci h, t + h)= n*(r, t) = m(r, t) + Ai (n(r, *)) . (34) 

The mass and momentum densities, p and j, are moments of the n^s with respect to 
the velocity vectors, 

p(r,t) = ^n i (r,i), (35) 

i 

j(r,*) = ^rii(r,t)ci, (36) 

i 

and therefore, the collision operator must satisfy the constraints of mass and momen- 
tum conservation, 

^A = ^Aci=0. (37) 

i i 

The LB algorithm has both locality and conservation laws built in, but two impor- 
tant symmetries have been lost. The system will in general exhibit cubic anisotropy, 
due to the underlying lattice symmetries, and violate Galilean invariance, due to the 
finite number of velocities. Isotropy can be restored in the large-scale limit by a 
careful choice of velocities and collision operator; however, the broken Galilean in- 
variance restricts the method to flows with u < c,. Since the speed of sound c s , 
the maximum velocity with which any signal can travel through the system, is of the 
order of the a, the condition actually means low Mach number (Ma) flow, 

Ma = u/c s < 1. (38) 

In soft-matter applications, variations in fluid density are small and there is a uni- 
versal equation of state characterized by the pressure at the mean fluid density and 

1 /2 

temperature, p = p{po, T), and the speed of sound c s — (dp/ dp) ' [42], 

p = Po + (p- p )c 2 s . (39) 

Within an unimportant constant (p — Pqc 2 s ), Eq. 39 can be replaced by the relation 
P = pel (40) 

which fits well to the linear structure of Eq. 34. The value of c s is immaterial except 
that it establishes a time-scale separation between sound propagation and viscous 
diffusion of momentum. For this reason, a model where c s is unphysically small 



16 



Burkhard Diinweg and Anthony J. C. Ladd 



may be used, so long as the dimensionless number C n = pcj/r] is sufficiently large; 
here I is a characteristic length in the system and r\ is the shear viscosity of the fluid. 
For polymers and colloids, C v <~ 10 — 1000, but values of C v in excess of 10 lead to 
quantitatively similar results. 

The simplest equation of state of the form of Eq. 40 is an ideal gas, 

p = ^k B T, (41) 

where T is the absolute temperature and k B Boltzmann's constant. Comparison with 
Eq. 40 yields 

k B T = m p c 2 8 . (42) 

The temperature is then determined by choosing values for the discretization param- 
eters b and h (c s ~ b/h), and the LB particle mass m p . The parameter m p controls 
the noise level in stochastic LB simulations [102]: the smaller m p (at fixed c s ), the 
smaller the temperature (or the noise level). This makes physical sense, since small 
m p means that a fixed amount of mass pb 3 is distributed onto many particles, and 
therefore the fluctuations are small. 

In this section we will study the connection between the LB equation, Eq. 34, 
and the equations of fluctuating hydrodynamics [42], 

dtp + d aJa = 0, (43) 
dtj a + dp (pc 2 s S al3 + pu a up) = df3<T al3 + dpa^p. (44) 

The Greek indexes denote Cartesian components, S a p is the Kronecker delta, and 
the Einstein summation convention is implied. The viscous stress has a Newtonian 
constitutive law, 

= Vafj-ySdyUg, (45) 

and for an isotropic fluid 

with shear and bulk viscosities r\ and rj v . The fluctuating stress tensor, a^ g , is a 



a/3' 

Gaussian random variable characterized by zero mean, { a^p ) = 0, and a covariance 



matrix 

(vie (r» *) <ts ( r '> *0) = 2k B Tr, a0l5 S (r - r') S (t - t') . (47) 



In the limit that T — > 0, cr^ vanishes, and the Navier-Stokes equations are recov- 
ered. 

We begin our analysis with a general description of the dynamics of the lattice- 
Boltzmann equation, based on a Chapman-Enskog expansion (Sec. 3.2). Then we 
consider the equilibrium distribution for the D3Q19 model (Sec. 3.3), followed by 
deterministic (Sec. 3.4) and stochastic (Sec. 3.5) collision operators. Finally we con- 
sider the connection of the fluctuating LB model to statistical mechanics (Sec. 3.6) 
and the effects of external forces (Sec. 3.7). 
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3.2 



Chapman-Enskog expansion 

The Navier-Stokes description of a fluid is more coarse-grained than the original 
LB equation, and to connect the microscopic scales with the hydrodynamic scales 
we follow a standard asymptotic analysis [103]. We first introduce a dimensionless 
scaling parameter e«l and write 



The idea is to measure spatial positions with a ruler that has such a coarse scale that 
details at the lattice level are not resolved. The position ri then corresponds to the 
number read off from this coarse-grained ruler; for example instead of talking about 
1000 nanometers, we talk about one micrometer. For two points to be distant on the 
hydrodynamic scale, it is not sufficient that \Ar\ is large, but rather that |Z\**i| is 
large. However, from the perspective of practical computation, the degree of coarse 
graining is never as extensive as implied by our analysis; the calculations would 
take far too long. Instead there is usually only a few grid points separating the lattice 
scale from the smallest hydrodynamic scale. Surprisingly the LB method can be quite 
accurate, even in these circumstances [99, 104]. 

In a similar way, we can also introduce a coarse-grained clock for the time vari- 
able, and write 



The fact that we choose the same factor e for both space and time is related to the 
typical scaling of wave-like phenomena, where the time scale of a process is lin- 
early proportional to the corresponding length scale. However, hydrodynamics also 
includes diffusion of momentum, where the time scale is proportional to the square 
of the length scale. These processes occur on a much longer time scale, and to capture 
the slow dynamics we introduce a second clock that is even more coarse-grained, 



We can therefore distinguish between "short times" on the hydrodynamic scale, char- 
acterized by t s = t\/e, and "long times", where ti — t 2 /e 2 . Both t s and ti are 
implicitly large on the lattice scale, with the hydrodynamic limit being reached as 
e — > 0. But once again, practical computation limits the separation between the time 
scales h, t s , and t\ to one or two orders of magnitude each. 

In the "multi-time scale" analysis, the LB population densities may be con- 
sidered to be functions of the coarse-grained position and times, ri, t\, and t 2 \ 
n.i = rii (ri,ti,t 2 )- When the algorithm proceeds by one time step, t — > t + h, 
ti — » <i + eh, and t 2 — > t 2 + e 2 h. The LB equation in terms of the coarse-grained 
variables is then, 

n,(r 1 + ec { h, h + eh, t 2 + e 2 h) - n;(ri, ii, i 2 ) = A, (n(ri, t u t 2 )) ■ (51) 

The population densities are slowly varying functions of coarse-grained vari- 
ables, and we may obtain hydrodynamic behavior by a Taylor expansion of m (Eq. 
51) to second order in powers of e: 



ri = sr. 



(48) 



t\ = et. 



(49) 



t 2 = e 2 t. 



(50) 
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n t (x + <5x) = n, (x) + ^ g^ Sxk + 2 ^ dx k dxf Xk ^ Xl + " ' ' ^ 

where we use x to indicate the coarse-grained variables, [ri, t\, t2\. Since the distri- 
bution function itself depends on the degree of coarse-graining, we must take the e 
dependence of the rii and Ai into account as well: 

n t =nf ) +en { p + 0(e 2 ), (53) 

Ai = 4 0) +eA { p +e 2 Af ] +0(e 3 ). (54) 

The conservation laws for mass and momentum must hold independently of the value 
of e, and thus at every order k: 

£A W =£4 fe) c<=0. (55) 

i i 

Inserting these expansions into Eq. 51, and collecting terms at different orders of e, 
we obtain: 

• at order e°, 

A (0) = 0; (56) 

• at orders 1 , 

(d tl +c 4 .d ri K (0) = fe- 1 A (1) ; (57) 

at order e 2 , 

dt 2 nf ] + ± (d tl + a ■ d ri ) 2 nf + (d tl + Cl ■ d ri )n^ = ft" 1 A {2) . (58) 

Subsequently, it will prove useful to eliminate the second occurrence of from 
Eq. 58, by using Eq. 57: 

3 t2 n^ + \ (d tl + c t ■ d Fl ) + n«) = h^Af \ (59) 

where n* = rn + Ai is the post-collision population in direction i. 

The multi-time-scale expansion of Eq. 5 1 is based on the physical time-scale 
separation between collisions it ~ h), sound propagation (t ~ h/e), and momentum 
diffusion (t <~ h/e 2 ). Eqs. 56-58 make the implicit assumption that these three relax- 
ations can be considered separately, which allows the collision operator at order k + 1 
to be calculated from the distribution functions at order k. In essence, the collision 
dynamics at order k + 1 is slaved to the lower-order distributions. The zeroth-order 
collision operator must be a function of n(°) only, 

A^ = A(n(°>), (60) 

which, in conjunction with Eq. 56, shows that n^ ) is a collisional invariant; thus 
we can associate n( ) with the equilibrium distribution n eq [105]. In order to avoid 
spurious conserved quantities, the equilibrium distribution should be a function of 
local values of the conserved variables, p and j, only. In a homogeneous system, with 
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fixed mass and momentum densities, n e9 (p,j) = n^(p, j) is stationary in time. A 
stochastic collision operator (see Sec. 3.5) cannot satisfy Eq. 56 and therefore must 
enter the Chapman-Enskog expansion at order e. 

From Eq. 53 we can derive analogous e expansions for p and j, 

p = p^+ep^+e 2 pW+0(e 3 ), (61) 

j=j( >+ £ j«+e 2 j( 2 )+O(e 3 ). (62) 

However, inserting these expansions into nj " 1 (p, j), shows that 

= = p^ =..., (63) 

0=j (1) =j (2) = ...; (64) 

otherwise nf" 1 would have contributions of order e and above, in contradiction to Eq. 
56. The mass and momentum densities can therefore be defined as moments of the 
equilibrium distribution as well, 

£< 9 = P, (65) 

i 

£«r c »=j- ( 66 ) 

i 

We can analyze the dynamics of the LB model on large length and time scales by 
taking moments of Eqs. 57 and 59 with respect to the LB velocity set Cj. From the 
zeroth moment, J2i ' " "> we obtain the continuity equation on the t\ time scale (Eq. 
55), 

d tl p + d la j a = 0, (67) 

and incompressibility on the ti time scale (Eqs. 55, 63, and 64) 

d t2 p = 0. (68) 

In Eq. 67 we have used the shorthand notation d\ a for the a component of the spatial 
derivative d Tl ■ 

The first moment, ^\ ■ ■ ■ Ci a , leads to momentum conservation equations on both 
time scales (Eqs. 55, 63, and 64): 

d t J a +d w ^l =0, (69) 
dt^ + ^fj^ =0, (70) 

where 7r Q( g is the momentum flux or second moment, 1 

TTgff = ^njCiqCip. (71) 

i 

Momentum is conserved on both the t\ and ti time scales, because, in the hydro- 
dynamic limit, the coupling between acoustic and diffusive modes is very weak. 



1 There is a notational inconsistency in Ref. [102]. In Eqs. 71 and 73 of that paper the super- 
script neq should be replaced by a superscript 1, and in Eq. 79 Q a p should be Q^. 



20 



Burkhard Diinweg and Anthony J. C. Ladd 



First, sound waves propagate with negligible viscous damping; then the residual 
pressure field in a nearly incompressible fluid relaxes by momentum diffusion. We 
can write the conservation laws on each time scale separately, as in Eqs. 69 and 
70, or combine them into a single equation in the lattice-scale variables r = ri/e 
and t = t\/e = t 2 /e 2 . The hydrodynamic fields depend on r and t parametrically, 
through their dependence on the coarse-grained variables n, t\, t%. Using d a for a 
component of d r , we have 

d a =edi a , (72) 

dt =ed tl +e 2 d t2 . (73) 

The combined equations for the mass and momentum densities on the lattice space 
and time scales are then: 

d t p + d a j a = 0, (74) 

dtja + drfb + \df> (C/ 9 + <7) = 0. C75) 

where from Eq. 53, tt^ = and tt™^ 9 = stt^ . 

Finally, we can derive a relation between the pre-collision and post-collision mo- 
mentum fluxes, 7r Q( g and -K* a(i , by taking the second moment of Eq. 57: 

du^+0^ = h^(^-^), (76) 
where <? Q( g 7 is the third moment of the distribution, 

$a(3- t = ^ n i c iaCi0C l7 . (77) 
i 

We note that is a collisional invariant and therefore remains unchanged by the 
collision process. In terms of the lattice variables, 

<T = <7 + h (d t nZ + a.**) . (78) 



'a/3 —"a/3 ' ^ * a/3 7 a/?7 

Equation 74 shows that continuity (Eq. 43) is automatically satisfied by any LB 
model. The Navier-Stokes equation (Eq. 44) will be satisfied, ;/ we succeed in en- 
suring that the Euler stress pc 2 8 a /3 + pu a up, the Newtonian viscous stress, a a p (Eq. 
45), and the fluctuating stress (Eq. 47) are given correctly by the sum of the 
momentum fluxes in Eq. 75. Since ir^ depends only on p and j, it must be identified 
with the Euler stress: 

^3 = P C l 5 a/3 + PUaUp. (79) 

The viscous stress and fluctuating stresses must then be contained in (?r*^ e9 + 
<?)/2- 

This is about as far as we can go in complete generality. In order to proceed fur- 
ther we need to consider specific equilibrium distributions and collision operators. 
The results of this subsection suggest the following approach towards constructing 
an LB method which (asymptotically) simulates the fluctuating Navier-Stokes equa- 
tions: 



Lattice Boltzmann simulations of soft matter systems 



21 



Find a set of equilibrium populations n e j q such that 

]Tnf = p; (80) 

i 

5>fc 4 =j; (81) 

i 

^2 n t q CiaCi(3 = pc 2 s S a p + pu a up. (82) 

i 

Find a collision operator Z\, with the properties: 

]TZ\, = 0; (83) 

i 

Ac, = 0; (84) 

i 

- The nonequilibrium momentum flux (tt*^ 69 + 7r™^' 3 )/2 must be connected 
with the sum of viscous and fluctuating stresses. 

In the following subsections, we will follow this procedure for the three-dimensional 
D3Q19 model. 

3.3 

D3Q19 model I: Equilibrium populations 

Early LB models [78-80] inherited their equilibrium distributions from LGCA, along 
with macroscopic manifestations of the broken Galilean invariance: an incorrect ad- 
vection velocity and a velocity dependent pressure. Subsequently a new equilibrium 
distribution was proposed that restored Galilean invariance at the macroscopic level 
[97, 106], but with the loss of the connection to statistical mechanics. The idea was to 
ensure that the first few moments of n q matched those derived from the Maxwell- 
Boltzmann distribution for a dilute gas [105], 

specifically; 

d 3 cn(c) = p, (86) 
d 3 cn(c)c a = pu a , (87) 



I- 

!■ 
I 



d 3 cn(c)c a Cf3 = — —S a p + pu a Uji = pc 2 s S a /3 + pu a up. (88) 

With these moments the Euler hydrodynamic equations (c.f. Eqs. 43 and 44), 

dtp + d a j a = 0, (89) 

dtja + dp (pc 2 s 5 a [i + pUaUjj) = 0, (90) 

may be derived from the continuum version of the Chapman-Enskog expansion 
[105]. The viscous stress arises from the non-equilibrium distribution (c.f. Sec. 3.2). 
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An expansion of the Maxwell-Boltzmann equilibrium distribution (Eq. 85) at 
low velocities suggests the following ansatz [97, 106] for the discrete velocity equi- 
librium, 

nf (p, u) = a c 'p (1 + Au- c t + B(u • c,) 2 + Cu 2 ) (91) 

with suitably adjusted coefficients a Ci , A, B, and C. The rationale for Eq. 91 is 
that the equilibrium momentum flux is quadratic in the flow velocity u (Eq. 88); it 
therefore makes sense to construct a similar form for n\ q . A drawback of Eq. 91 
is that n eq may be negative if u becomes sufficiently large. This can be avoided by 
more general equilibrium distributions, which are equivalent to Eq. 91 up to order u 2 
[107-109]. 

The prefactors a c% > are normalized such that 



(92) 



which ensures that Eq. 65 is satisfied in the special case u = 0. The notation in Eq. 
91 was chosen in order to indicate explicitly that the weights depend only on the 
absolute value of the speed c», but not its direction; this follows from the rotational 
symmetries of the LB model. The coefficients A, B, C are here independent of Cj. 
There are other LB models, like D3Q18 [106], where this condition is not imposed, 
and A, B, and C depend on as well; however, such models are only hydrody- 
namically correct in the incompressible limit [98], and cannot be straightforwardly 
interpreted in terms of statistical mechanics (see Sec. 3.6). We will not consider such 
models. 

In a cubic lattice, symmetry dictates the following relations for the low-order 
velocity moments of the weights, 

5> Cl c M =0, (93) 

i 

a Ci CiaCi/3 = C 2 5 af3 , (94) 

i 

^2a c *c la c lf 3C l7 = 0, (95) 

i 

^ a Cz c ia Cif3C i7 Cis = C' A S al 3js + C4 (S a/3 8 7 s + 5 a7 6ps + dasdpy) , (96) 

i 

where the values of the parameters C 2 , C4 and C 4 depend on the details of the choice 
of the coefficients a Ci . The tensor S a ^s is unity when a = j3 = 7 = S and zero 
otherwise. Rotational invariance of the stress tensor requires that C' 4 = 0. 

The results in Eqs. 92-96 allow us to calculate the moments of Eq. 91 up to 
second order. Consistency with the mass density, momentum density and Euler stress 
for a given p and u, uniquely determines the equilibrium distribution, 
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with the speed of sound c 2 = C 2 , and the weights adjusted such that C4 = and 
C4 = Cf- These two latter conditions, together with the normalization condition 
Eq. 92, form a set of three equations for the coefficients a Ci . Therefore at least three 
speeds, or three shells of neighbors, are needed to satisfy the constraints. We consider 
the D3Q19 model, which incorporates the three smallest speeds on a simple-cubic 
lattice. Here one obtains a = 1/3 for the stationary particles, a 1 = 1/18 for the six 
nearest-neighbor directions, and a^ 2 = 1/36 for the twelve next-nearest neighbors. 
The speed of sound is then c 2 s = (1/3) (b/h) 2 . 

We now turn back to the results of the previous subsection, since the explicit form 
of n^ q allows us to pursue the analysis further. We first calculate the equilibrium 
third-order moment (Eq. 77) using Eq. 97: 

®ap-y = P C l ("a Vr + u ^ al + W 7 <5 Q/3 ) . (98) 

In fact Eq. 98 is model independent to order u 2 , since only the linear term in u • Ci 
contributes. To close the hydrodynamic equations for the mass and momentum den- 
sities (Eqs. 74 and 75) we need expressions for the pre-collision and post-collision 
momentum fluxes, 7r*^ 69 and ■ From Eq. 76 we can obtain an expression for 
^*a^ q ~ m terms °f th e velocity gradient, 

pc 2 s (d la u + d w u a ) = h- 1 (tt*^ - tt$) , (99) 

where we have used Eqs. 67 and 69 to rewrite the time derivative of 7r^ in terms of 
spatial derivatives of p and u. In arriving at Eq. 99, we have neglected terms of order 
u 3 , consistent with the low Mach number limit we are considering. Finally, Eq. 99 
can be rewritten in terms of the unsealed variables, 

tt*J 9 - tt™^ = hpc 2 s (d a u p + dpu a ) . (100) 

To obtain a further relation for the non-equilibrium momentum fluxes, we must con- 
sider the collision operator in more detail. 



3.4 

D3Q19 model II: Deterministic collision operator 

In a deterministic model, the collision operator Ai is a unique function of the distri- 
bution n. Therefore, we can obtain the Chapman-Enskog ordering of Ai via a Taylor 
expansion with respect to n: 



Ai (n) = Ai (n 

(■<«>) + .£(f£ 



l(0) 



nf > + O (e 2 ) . 



(101) 



The analysis of Sec. 3.2 has shown that A(n (0) ) = (Eq. 60), and that hydrody- 
namic behavior is determined by the order e 1 collision operator, 
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(i) _ v / dAj 
i 4^ V dn 4 



nf\ (102) 



(2) 

Although A\ ' appears at second order in the Chapman-Enskog expansion (Eq. 58), 
it makes no contribution to the change in mass and momentum densities (Eq. 55). 
However, A^p contributes to a first-order change in the viscous stress (Eq. 78), 
which enters into the momentum equation at second order (Eq. 70). It is therefore 
reasonable to construct the collision operator with the form of A^ : 

A i = Y J £ijnY\ (103) 

3 

where dj is a matrix of constant coefficients. Thus to lowest order in e, the collision 
process is a linear transformation between the non-equilibrium distributions for each 
velocity: 

«r=E^ +£ «) n r- (io4) 

3 

The simplest such collision operator is the lattice BGK (Bhatnagar-Gross-Krook) 
model [77], Cij = —Sij/r, where the collisional relaxation time r is related to 
the viscosity. Here we will work within the more general framework of the multi- 
relaxation time (MRT) model [110], for which the lattice BGK model is a special 
case. 

Polynomials in the dimensionless velocity vectors, c, = Cj/c (c = b/h), form 
a basis for a diagonal representation of Cij [110], which allows for a more general 
and stable LB model with the same level of computational complexity as the BGK 
version [111]. Orthogonal basis vectors, e^, are constructed from outer products of 
the vectors Cj. For example: 

e * = 1, (105) 

eu = c ix , (106) 

e2i = c iy , (107) 

e 3l =c lz . (108) 

There are six quadratic polynomials, which are given in Table 1 as basis vectors 
e 4 — e g . A Gram-Schmidt procedure ensures that all the basis vectors are mutually 
orthogonal with respect to a set of positive weights, q Ci > 0, 

^2 1 Cze kteu = w k Ski- (109) 

i 

The weights are restricted by the same symmetries as the coefficients in the equilib- 
rium distribution a Ci , but are not necessarily the same; in the D3Q19 model there are 
then three independent values of q Ci . The normalization factors, Wk > 0, are related 
to the choice of basis vectors 

w k = J2^ e l- (HO) 
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Table 1 . Basis vectors of the D3Q19 model. Each row corresponds to a different basis vector, 
with the actual polynomial in c ia = c ia /c shown in the second column. The normalizing 
factor for each basis vector is in the third column. The polynomials form an orthogonal set 
when q c ' = a c ' (Eq. 109). 
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Cki 







l 


1 


1 




1/3 
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Ciy 


1/3 


3 




1/3 
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c? - 1 


2/3 


5 




4/3 


6 


r 2 r 2 

Ciy Wz 


4/9 


7 


Cix Ciy 


1/9 


8 


CiyCiz 


1/9 


9 


CizCix 


1/9 


10 


(3c^ — 5)cj x 


2/3 


11 


(3c^ 5) 
(3c? - 5)c lz 


2/3 


12 


2/3 


13 




2/9 


14 




2/9 


15 




2/9 


16 


3c* - 66? + 1 


2 


17 


(2c?-3)(3cL-c?) 
(2cf-3)(4-cL) 


4/3 


18 


4/9 



Within the D3Q19 model, polynomials up to second order are complete, but at third- 
order there is some deflation; for example, cf x is equivalent to Ci X . In fact, there are 
only six independent third-order and three independent fourth-order polynomials in 
the D3Q19 model. Beyond fourth order, all polynomials deflate to lower orders, so 
the basis vectors in Table 1 form a complete set for the D3Q19 model. 

The basis vectors can be used to construct a complete set of moments of the LB 
distribution, 

77J fc = ^ e fc ;ri;, (111) 

i 

which allows for a diagonal representation of the collision operator [110, 112], as 
will be made clear later. Hydrodynamic variables are related to the moments up to 



quadratic order in £j (c.f. Table 1): 

p = m , (112) 

j x =m 1 c, (113) 

j y =m 2 c, (114) 

jz = m 3 c, (115) 

n xx = (mo + m 4 + m 5 )c 2 /3, (1 16) 
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^yy 


= (2m - 


h 2m,4 




= (2m H 


- 2m4 


T^xy 


= myc 2 , 




Kyz 


= m 8 c 2 , 




Kzx 


= mgc 2 . 





m 5 + 3m 6 )c 2 /6, (117) 
m 5 - 3to 6 )c 2 /6, (118) 

(119) 
(120) 
(121) 

There are additional degrees of freedom in the D3Q19 model beyond those required 
for the conserved variables and stresses (Eqs. 112-121). These "kinetic" or "ghost" 
[101] moments do not play a role in the large-scale dynamics [102], but they are 
important for proper thermalization [101] and near boundaries [113]. 

The basis vectors in Table 1 are complete but not unique. Besides trivial varia- 
tions in the Gram-Schmidt orthogonalization, there is a substantive difference that 
depends on the choice of the weighting factors q Ci : these factors determine both the 
result of the orthogonalization procedure, as well as the back transformation from 
moments m k to populations rn. This is most easily seen from the observation that 
Eq. 109 can be rewritten as the standard orthonormality relation [102] 

y^efc»ea = 5 k i, (122) 

i 

where we have introduced the orthonormal basis vectors 

eki = \ — e ki . (123) 
V w k 

Equation 122 implies the backward relation 

ekiikj = Sij, (124) 

k 

or, in terms of unnormalized basis vectors, 

y^Wfc 1 efcje fc j = -\-5i - (125) 

k q * 

In the normalized basis the transformations between distribution and moments are 

m k = —= = } e ki —= = > e ki rii, (126) 

ni = —= = }ehi—= = y^e ki m k ; (127) 
Vq^ Y ^ k 

we will make use of these relations in Sec. 3.5. The analog of Eq. 104 for the nor- 
malized basis is 

«r e9 = E(^+4H ne9 ' d28) 

3 

with 
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In terms of unnormalized basis vectors the back transformation is given by 

m = q Ci ^2 w~ k 1 e kl m k . (130) 
k 

The most obvious choice is to set q Ci = 1 [110, 112], but then the basis vectors 
of the kinetic modes, eio — eis, are not orthogonal to the equilibrium distribution, 
and the moments m u — mis nave both equilibrium and non-equilibrium contri- 
butions [110, 112]. The statistical mechanical connection is more straightforward 
if the weights q Ci are matched to the weights in the equilibrium distribution, setting 
qCi _ a a . tn j s eliminates the projection of n 9 on the kinetic moments. The weighted 
orthogonality relation defines a different but equivalent set of basis vectors to those 
given in Refs. [110, 1 12], and these are the ones given in Table 1. A comparison of 
the two sets of basis vectors can be found in Ref. [1 14]. 

The basis vectors can be used to construct a collision operator that automatically 
satisfies all the lattice symmetries, 

A? = ^ Afcgfciefcj, (131) 

k 

which is a symmetric matrix, while Cij, in general, is not symmetric. The orthogo- 
nality of the basis vectors ensures that each moment relaxes independently under the 
action of the linearized collision operator, 

rnl ne i = lk m n k e \ (132) 

where j k — 1+A&. For the conserved modes k = 0, . . . , 3 the value j k is immaterial, 
since m^ 9 = rn^ eq = 0. For the other modes, k > 3, linear stability requires that 

hk\ < i; (133) 

i.e. the effect of collisions must be to cause the nonequilibrium distribution to de- 
crease rather than increase. The eigenvalues, 7^, may be positive or negative, with 
7/t < corresponding to "over-relaxation". 

The number of independent eigenvalues is limited by symmetry. There are at 
most six independent -f^s in the D3Q19 model, corresponding to a bulk viscous 
mode with eigenvalue j v , five symmetry-related shear modes, which must have the 
same eigenvalue, 7 S , and nine kinetic modes, broken down into symmetry-related 
groups: eio — ei 2 , ei3 — ei5, ei 6 , en — ei 8 . The eigenvalues -f s and -y v can be related 
to the shear and bulk viscosities by decomposing the stress tensor into traceless- 
symmetric (shear) and trace (bulk) components, 

= K aj 3 + -7r 77 <5 Q/ 3; (134) 

the overbar is used to denote a traceless tensor. Equation 132 implies the following 
relations between pre- and post-collisional stresses: 

T *neq = ^ T ne 9 _ (136) 
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Additional relations between the pre- and post-collision stresses have already been 
provided (Eq. 100): 

<7 Q ~ K7 = hpcl (9^ + d^ a ) , (137) 

KT ~ <l q = 1hpc 2 s d a u a . (138) 

Equations 135-138 can be solved to relate the pre- and post-collision stresses to the 
velocity gradient: 

_ hP± { ^ + Q^ ah (139) 

1 Is 
1 Is 

1 — 7f 

< n a eq = -^^d a u a . (142) 

From Eq. 75 we then find the usual Newtonian form for the viscous stress, and can 
identify the shear and bulk viscosities: 

Lattice symmetry dictates that there are at most four independent eigenvalues of 
the kinetic modes: (see Table 1): 7 3a (modes 10— 12), 73;, (modes 13—15), 74a (mode 
16), and 746 (modes 17 — 18). In a number of implementations of the MRT model 
[99, 100, 106, 115] the kinetic eigenvalues are set to zero, so that these modes are 
projected out by the collision operator, although they reoccur at the next time step. 
Recently, it has been shown that the kinetic eigenvalues can be tuned to improve the 
accuracy of the boundary conditions at solid surfaces [113]. A useful simplification 
is to use only two independent relaxation rates, with j v = -f s = j ia = 74b = j e 
and 7 3a = 73b = j a . The optimal boundary conditions are obtained with specific 
relations between -f e and j [113, 1 14]. 



3.5 

D3Q19 model III: Thermal noise 

In the fluctuating LB model [98, 100], thermal noise is included by adding a stochas- 
tic contribution, A' i7 to the collision operator: 

The collision operator must still conserve mass and momentum exactly, 
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E4 = E4« 



(146) 



while the statistical properties of A[ include a vanishing mean, (Z\-) = 0, and a 
nontrivial co variance matrix, (Z\^Z\^), that gives the correct fluctuations at the hy- 
drodynamic level (see Eqs. 44 and 47): 



2k B T 
b 3 h 



-Va/3-rS- 



(147) 



The stochastic collision operator is assumed to be local in space and time, so that 
there are no correlations between the noise at different lattice sites or at different 
times. The delta functions in space and time have been replaced by b~ 3 and h~ x , 
respectively, so that the double integral of Eq. 47 with respect to r' and t', over a 
small space-time region of size b 3 h, matches the corresponding integral of Eq. 147. 

Splitting the tensor into the trace, a^ a , and traceless, a f aj3 , parts gives the equiv- 
alent relations 



a a f3 a j6 



if 

T af3 l 



rf 



2k B T V 
b 3 h 



b 3 h 



= 0. 



(148) 

(149) 
(150) 



Although temperature does not appear directly in the D3Q19 LB model, we can de- 
termine the appropriate fluctuation level through the equation of state for an isother- 
mal ideal gas of particles of mass m p , k B T = m p c 2 s — (J>b 3 c 2 [102]. Taking into 
account the results for rj and r\ v (Eqs. 143 and 144), we can write the desired corre- 
lations in terms of the LB variables: 



-/ -/ 

T ap (J 1 8 



l + 7.s 
1-7, 

6- 1 



l-7t> 



0. 



(151) 

(152) 
(153) 



The stress fluctuations a^p 316 different from the random stresses o r aS that arise 



in the LB algorithm itself, 



a/3 



'a/3 



(154) 



The reason is that cr^ pertains to fluctuations on the t\ time scale, which interact 
with the hydrodynamic flow field, while <j r a „ represents added noise on the lattice 
time scale, h. We use the Chapman-Enskog procedure to work backwards from the 
known fluctuations in a f aj3 to determine the covariance matrix for a r a p [102]. The 
stress update rule including random noise a r aj3 is (c.f. Eqs. 135 and 136): 
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*a~P = ~T=t ( 9aU P + + T^V^' (15?) 



-*neq - neq . - r CC x 

* a fi -Is* a? ( 155 > 

<o e9 = 7»0 + <£«- (156) 

Eqs. 137 and 138 remain valid and, together with Eqs. 155 and 156, can be solved 

for the pre- and post-collisional stresses, 7r™^ 9 and 7r*^ e<? , as before: 

IS 1 

<T = -^r 1 + W + 7^-^ d58) 

< e a q = -^ L d a u a + -^—a r aa , (159) 
1 — 7t» 1 — 7« 

<neq = JhP±hL daUa + 1 a r (160) 

1 - 7w 1 - 7t> 

Comparing Eq. 75 with Eq. 44 we can read off the relations between the hydrody- 

namic fluctuations and the random noise, 

1 

'a/3 = ~JZ ~ U af3> 
Is 

°L = —, (162) 



^ aB = -T^-<>- (161) 
1 



1 - 7« 

Therefore, the random noise inserted at the microscopic (LB) level must have the 
following covariances: 



°a/3 c V 



^4 



(163) 



= 6(1- 7 ;) , (164) 

= 0. (165) 

The random stress has a typical amplitude of y/Jipc^ and is obtained from the 
second-order moment of the fluctuations in n" e9 . Therefore a typical fluctuation in 
the population density is of order y/JIp. Combining this scaling with Eqs. 126 and 
127, suggests dimensionless variables 

fit = — — , (166) 

m k = , (167) 

which transform using the symmetric basis vectors defined in Eq. 123, 

m k = ^e k ihi, (168) 

% 

nj = y^e k iifi k . (169) 

k 
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The stochastic collision operator can then be implemented independently for each 
mode m k , 

mT q =lkm n k eq + Vkr k , (170) 

where r k are independent Gaussian random variables with zero mean and unit vari- 
ance. The dimensionless constants if/, are determined by expressing the random 
stresses a T a ^ in terms of the r k and cp k and then calculating the covariance matrix 
[102]. For example, u r = ^/JIpc^(pir 7 , while a r aa — V6mP c sV 3 4 7 '4- Comparison 
with Eqs. 163 - 165 shows that the correct stress correlations are obtained for 

Vk = {l-lt) 1/2 . (171) 

At the hydrodynamic scale, only fluctuations in stress contribute to the time evo- 
lution of the momentum density (Eq. 44) so in principle it is sufficient to add random 
fluctuations to the modes m^, . . . ,mg only: In the original derivation of the fluc- 
tuating LB equation [98, 100], the kinetic modes were projected out entirely, i.e. 
Ik = fk = for k = 10, 11, ... , 18. More recently, Adhikari et al. [101] have 
argued that the kinetic modes should be thermalized as well. They extended Eq. 170 
to the kinetic modes (k = 10, . . . , 18), with -y k = (as in Refs. [98, 100]) but with 
fk = 1, which then satisfies Eq. 171. It was demonstrated numerically that this leads 
to more accurate fluctuations at short length scales, but the theoretical justification 
remained somewhat obscure. From the discussion so far, we can see that both pro- 
cedures give the same random stresses <J r a p, and hence are not different from the 
point of view of fluctuating hydrodynamics. This has been clarified recently [102], 
by analyzing the LB model in terms of statistical mechanics. A purely microscopic 
approach was taken, in which the stochastic collisions were viewed as a Monte Carlo 
[116] process. Knowledge of the probability distribution of the LB variables n then 
makes it possible to check whether or not a given collision rule satisfies the condition 
of detailed balance. It can be shown [102] that the kinetic modes must be thermalized 
in order to satisfy detailed balance, in agreement with the procedure proposed in Ref. 
[101]. The theory will be outlined in Sec. 3.6. 



3.6 

Statistical mechanics of lattice-Boltzmann models 

The starting point of the statistical mechanical development in Ref. [102] is the no- 
tion of a generalized lattice gas. We define fj(r, t) in Eq. 30 as the number of par- 
ticles with velocity c, at site r at time t. In contrast with the standard LB model, 
Vi is a (positive) integer; in contrast with lattice-gas models, i/j » 1. The state at 
a particular lattice site, i/(r, t), is modified by the collision process, subject to the 
constraints of mass and momentum conservation; the post-collision state, v*(r,t), 
is then propagated to the neighboring sites (Eq. 30). 

Although a deterministic GLG collision operator would be difficult to construct, 
we can nevertheless determine the distribution in a homogeneous equilibrium state 
from the conservation laws alone. First we note that there is an entropy associated 
with each v;, 
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Si = - (vi In Vi-Vi- Vi In v t + v{) , (172) 

where i>i is the mean value of Vi in the homogeneous state. Each velocity direction i 
at each lattice point has a degeneracy exp(S'j), which can be derived from a Bernoulli 
process. Particles are selected for the velocity direction i with probability po, with 
pa chosen so that on average a total of v = N r po particles will be selected from a 
reservoir of N r particles. Then the probability to select exactly v particles is given 
by the binomial distribution, 

N r l ( v\ v ( v \ Nr ~ v 

pM = ^v^UJ v-wj ■ (173) 

Equation 172 results from calculating lnp(^) in the limit of N r — > oo, at fixed v. 
Under the usual assumption that in the equilibrium state the populations correspond- 
ing to different lattice sites and different directions are uncorrected, the entropy per 
lattice site is Siy) = ^ Si. 

The populations at a given lattice site are sampled from a probability distribution 
proportional to exp [S (i/)], but subject to the constraints of fixed mass and momen- 
tum density, which characterize the homogeneous state: 

P (v) oc exp [S (v)] S^t^Ui-pjs^Yl v ^ - $J ■ ( 174 ) 

Consistency with the formalism developed in the previous sections requires 

fiDi — pa Ci . (175) 

The equilibrium or mean populations for a given p and j are found by maximizing 
P or, more conveniently, by maximizing S and taking into account the conservation 
laws by Lagrange multipliers: 

dS 

— + \ p + Aj • a = 0, (176) 

i 
i 

The exact solution is 

9 = Vi exp (Ap + Aj • Ci) , (179) 

where the Lagrange multipliers, A p and Aj, are found from the constraint equations 
177 and 178. Solving these equations in terms of a power series in u, and disregarding 
terms of order 0(u 3 ), one finds the standard equilibrium distribution given in Eq. 97. 
This approach has been previously proposed within the framework of the "entropic 
lattice-Boltzmann" method [108, 109], which however, focuses exclusively on the 
deterministic LB model. 

Within the statistical-mechanical framework we have developed for the LB 
model, the population densities rii fluctuate around mean values determined by the 
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hydrodynamic flow fields. Thus the non-equilibrium distribution is sampled from 
P(n neq ) which is Gaussian distributed about the equilibrium [102], 



p („-) oc exp ( - y: ) 6 ( £ < eq ) 6 ( E c < »r ) • < 18 °) 



The variance of the fluctuations is controlled by the mass density /j,, associated with 
an LB particle. A small number of particles gives rise to large fluctuations and vice 
versa. For simplicity we will ignore the effects of flow on the variance of the distribu- 
tion, replacing n eq by its u = value. This can be justified at the macroscopic level 
by the Chapman-Enskog expansion [102]. Rewriting Eq. 180 in terms of normalized 
variables hi (see Eq. 166), and transforming to the normalized modes (see Eqs. 167 
and 169) eliminates the explicit constraints, 

P(rh neq ) a exp ^gm^ 2 ) ■ (181) 

Fluctuations only arise in the non-conserved modes, while the conserved modes have 
no non-equilibrium contribution, i.e. m r k ieq — for k < 3. 
We now reinterpret the update rule Eq. 170, 

; *net/ ~ neq . /i O o\ 

m k = "/km k + (p k r k , (182) 

as a Monte Carlo move. The transition probability is then identical to the probability 
of generating the random variable r k , 

( - *neq ~ neq \ 2 



K e, -^n=(2V fc ) _1/ exp 



(183) 



and for the reverse transition the same formula holds, with the pre- and post- 
collisional populations exchanged. The condition of detailed balance [116], 

u (m n k eq -» m k neq ) _ exp (-(m k neq ) 2 /2) 



lu (m* k neq - m n k eq ) exp (-(m«^) 2 /2) ' 
then holds if and only if 

^ = (l-7fe) 1/2 ' (185) 

as before (Eq. 171). The important point is that this relation, which in the previ- 
ous subsection was only proved for the stress modes, can now be shown to hold for 
all non-conserved modes. It is a necessary condition for consistent sampling of the 
thermal fluctuations, not just on the macroscopic hydrodynamic level (for which the 
stress modes alone are sufficient), but also on the microscopic LB level itself. Al- 
though assigning j k = (and tpk = 1) to all kinetic modes is obvious and straight- 
forward [101], the present analysis shows that this is not necessary. Other values of 
7^ and ipk are possible as well, so long as they satisfy Eq. 185, and specific values 
may be desirable for a more accurate treatment of boundary conditions [113, 1 14]. 
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3.7 

External forces 

An external force density f (r, t) can be introduced into the LB algorithm by an ad- 
ditional collision operator A", 

3 

For simplicity, we only consider the deterministic case; the analysis of the fluctuating 
part (Zi-) remains the same. Application of the new collision operator should leave 
the mass density unchanged, but increase the momentum density by hi. This implies 
the following conditions on the moments of A"\ 

J2A? = 0, (187) 

i 

Y,A"ci = hf. (188) 

i 

Consequently the definition of the fluid velocity is no longer unique: one can legiti- 
mately choose any value for u between p^ 1 riiCi) and p^ 1 n i°i + hf) 0- e - 
between the pre-collisional and post-collisional states). However, numerical [98] 
and theoretical [100, 117, 118] analysis shows that the optimum value is just the 
arithmetic mean of the pre- and post-collisional velocities. We define the momen- 
tum density as 

j = ^n t c, + ^f, (189) 

i 

and the corresponding flow velocity as u = ]/ p. Consistency with Eq. 66 requires 
that we use this value for u to calculate n^ 9 (Eq. 97): 



]r< 9 c J= j, d90) 

i 

X)nr«c i = -^f. (191) 

i 

In Ref. [100] the usual moment condition J^. n™ eq Ci = was maintained. In com- 
parison with the present approach this makes a small error of order / 2 to the distri- 
bution, which leads to spurious terms in the Chapman-Enskog analysis. In contrast 
the present approach leads to a clean result, entirely equivalent to the force-free case. 
This may be of consequence when there are strongly inhomogeneous forces, such 
as are considered in Sec. 4. However, it should also be noted that these differences 
vanish in the low Reynolds number limit. 

Since A"(n eq ) ^ 0, the Chapman-Enskog expansion of A" starts at order e 1 
(c/Eq.56): 

A'! = eA'l {1) +e 2 A'l {2) + .... (192) 
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Following the procedure of Sec. 3.2, we take moments of Eqs. 57 and 59 and obtain 
similar equations for the mass and momentum density (c.f. Eqs. 74 and 75): 



dtp + d a j a = 0, 
d t j a + dpn^p + 



\df> 



*ncq i neq 
7T„, a + 7T 



'a/3 



'a/3 



/a- 



(193) 
(194) 



However, the second moment leads to a force-dependent contribution to the non- 
equilibrium momentum flux, which can be derived as before, beginning with Eq. 76 
and substituting the equilibrium expressions for (Eq. 79) and ^>^g 7 (Eq. 98). 
The time derivative of the momentum flux now generates terms involving uf + f u 
from the momentum conservation equation on the t\ time scale (c.f. Eq. 99): 

<J q - <7 = h Pcl (9 a u + d fj u a ) + h (u a f (i + up f a ) . (195) 

Spurious terms proportional to uf can be eliminated from Eq. 194, by including 
the second moment of A", 



(196) 



so that the stress update is now (c.f. Eq. 135 and 136) 



*neq 



1 



'a/3 = IsK? + pvTr™ Q S al 3 + hS af3 . 



(197) 



Equations 195 and 197 form a linear system for 7r™^ 9 and 7r*^ e<? . Solving these equa- 
tions as before (Eqs. 139 - 142), and inserting the result into Eq. 194, we obtain a 
Newtonian stress with unchanged values for the viscosities by choosing £ a/ 3 such 
that 



U a f(3 + Upf a - -M 7 / 7 (5a/3 



+ ^(l + 7v)«7M*/3- ( 198 > 



The moment conditions expressed by Eqs. 187, 188 and 198 are uniquely satisfied 
by the choice 



A 1 - 



h 



h 

2rf 



S aj 3 (ciaCi/3 - C 2 s 5 a p) 



(199) 



where A" only affects the modes mi, . . . , mg. This result has been derived pre- 
viously [118] within the context of the LBGK model; here we have presented the 
derivation in the more general MRT framework. 



4 

Coupling the LB fluid to soft matter 

The fundamental algorithmic problem in soft matter simulations is the coupling be- 
tween the solid and fluid phases. A key attraction of LB methods is the simplicity 
with which geometrically complex boundaries can be incorporated. The first correct 
implementation of a moving boundary condition was described in the proceedings 
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of a workshop on lattice-gas cellular automata [119]; a more accessible source is 
Ref. [120]. The idea was to modify the bounce-back rule for stationary surfaces such 
that the steady-state distribution was consistent with the local surface velocity. By 
constructing the boundary-node interactions along the individual links, the viscous 
stress remains unchanged. Subsequently, we showed numerically that this algorithm 
gives accurate hydrodynamic interactions between spherical particles suspended in 
a lattice-gas fluid [120]. Nevertheless, it quickly became clear that the fluctuating 
LB model was a more useful computational tool, for the reasons outlined in Sec. 3. 
The LG algorithm for the moving boundary condition carries over in a simple and 
direct way to the LB method [98]. A number of improvements to the bounce-back 
boundary condition have been proposed over the years, and we will summarize some 
of the more practical approaches in Sec. 4.4. 

More recently an entirely different approach has been proposed [121, 122] in 
which particles couple to the fluid through a frictional drag. This method has the 
advantage of greatly reducing the number of LB grid points in the simulation, at the 
cost of a representation that is only correct in the far field. The method has been ap- 
plied to polymers [43, 122, 123] and to suspended solid particles [124-128]. In the 
latter case the surface is described by a number of sources distributed over the surface 
of the particle. The distributed forces resemble the Immersed Boundary (IB) meth- 
ods [129], which are common in finite-difference and finite-element simulations; this 
connection has only been recognized recently [130]. We will summarize these devel- 
opments and add some new ideas and interpretation of the force coupling methods. 
In a related development [131, 132], conventional immersed boundary methods are 
being used in conjunction with an LB fluid. However, the coupling in this case is 
implicit, solving for the velocity of the interface through a force balance, which cor- 
responds to the high friction limit of Refs. [124-128]. Here, we will only consider 
inertial coupling, since the theory for thermal fluctuations has not been worked out 
for the implicit schemes. 

4.1 

Boundary conditions 

To simulate the hydrodynamic interactions between solid particles in suspension, the 
lattice-Boltzmann model must be modified to incorporate the boundary conditions 
imposed on the fluid by the solid particles. The basic methodology is illustrated in 
Fig. 1 . The solid particles are defined by a boundary surface, which can be of any size 
or shape; in Fig. 1 it is a circle. When placed on the lattice, the boundary surface cuts 
some of the links between lattice nodes. The fluid particles moving along these links 
interact with the solid surface at boundary nodes placed halfway along the links. 
Thus a discrete representation of the particle surface is obtained, which becomes 
more and more precise as the particle gets larger. 

In early work, the lattice nodes on either side of the boundary surface were treated 
in an identical fashion [98, 120], so that fluid filled the whole volume of space, both 
inside and outside the solid particles. Although the fluid motion inside the particle 
closely follows that of a rigid solid body [99], at short times the inertial lag of the 
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Fig. 1. Location of boundary nodes for a curved surface. The velocities along links cutting 
the boundary surface are indicated by arrows. The locations of the boundary nodes are shown 
by solid squares, and the fluid nodes by solid circles. The open circles indicate nodes in the 
solid adjacent to fluid nodes. 

fluid is noticeable, and the contribution of the interior fluid to the particle force and 
torque reduces the stability of the particle velocity update. Today, most simulations 
exclude interior fluid, although the implementation is more difficult when the parti- 
cles move. The moving boundary condition [98] without interior fluid [133] is then 
implemented as follows. We take the set of fluid nodes r just outside the particle sur- 
face, and for each node all the velocities c& such that r + Cf,/i lies inside the particle 
surface. An example of a set of boundary node velocities is shown by the arrows in 
Fig. 1 . Each of the corresponding population densities is then updated according to 
a simple rule which takes into account the motion of the particle surface [98]; 



where n\ (r, t) is the post-collision distribution at (r, t) in the direction and cy = 
—Cb. The local velocity of the particle surface, 



is determined by the particle velocity U, angular velocity fi, and center of mass R; 
rb = r + jhcb is the location of the boundary node. 

As a result of the boundary node updates, momentum is exchanged locally be- 
tween the fluid and the solid particle, but the combined momentum of solid and fluid 
is conserved. The forces exerted at the boundary nodes can be calculated from the 
momentum transferred in Eq. 200, and the particle forces and torques are then ob- 
tained by summing over all the boundary nodes associated with a particular particle. 
It can be shown analytically that the force on a planar wall in a linear shear flow is 
exact [98], and several numerical examples of lattice-Boltzmann simulations of hy- 
drodynamic interactions are given in Ref. [99]. Figure 2 illustrates the accuracy that 
can be achieved with the MRT collision operator described in Sec. 3.4. Even with 
small particles, only 56 in diameter, the hydrodynamic interactions are within 1% of 
a precise numerical solution [21], down to separations between the particle surfaces 
s = r — 2a ~ b, corresponding to s <~ 0.4a, where a is the sphere radius. Periodic 



nv(r,t + h) = n* b (r,t) 



2a Cb p\ib ■ Cb 



(200) 



u 6 = U + fix(r b -R), 



(201) 



38 



Burkhard Diinweg and Anthony J. C. Ladd 




Fig. 2. Hydrodynamic interactions from LB simulations with particles of radius a = 2.56. 
The solid symbols are the LB friction coefficients, (^ pU and (P er P^ for the relative motion of 
two spheres along the line of centers (left) and perpendicular to the line of centers (right). 
Results are compared with essentially exact results from a multipole code [21] in the same 
geometry (solid lines). 



boundaries with a unit cell size L — 12a were used, with the pair inclined at 30° to 
a symmetry axis; other geometries give a very similar level of agreement. We em- 
phasize that there are no adjustable parameters in these comparisons. In particular, 
in contrast to previous work [99, 104], there is no need to calibrate the particle ra- 
dius; the correct particle size arises automatically when the eigenvalues of the kinetic 
modes of the MRT model have the appropriate dependence on the shear viscosity 
[114]. 

To understand the physics of the moving boundary condition, one can imagine 
an ensemble of particles, moving at constant speed Cb, impinging on a massive wall 
oriented perpendicular to the particle motion. The wall itself is moving with velocity 
Ub <C Cft. The velocity of the particles after collision with the wall is — Cb + 2u& 
and the force exerted on the wall is proportional to — U;,. Since the velocities in 
the lattice-Boltzmann model are discrete, the desired boundary condition cannot be 
implemented directly, but we can instead modify the density of returning particles 
so that the momentum transferred to the wall is the same as in the continuous veloc- 
ity case. It can be seen that this implementation of the no-slip boundary condition 
leads to a small mass transfer across a moving solid-fluid interface. This is physi- 
cally correct and arises from the discrete motion of the solid surface. Thus during a 
time step h the fluid is flowing continuously, while the solid particle is fixed in space. 
If the fluid cannot flow across the surface there will be large artificial pressure gra- 



Lattice Boltzmann simulations of soft matter systems 



39 



dients, arising from the compression and expansion of fluid near the surface. For a 
uniformly moving particle, it is straightforward to show that the mass transfer across 
the surface in a time step h (Eq. 200) is exactly recovered when the particle moves 
to its new position. For example, each fluid node adjacent to a planar wall has 5 links 
intersecting the wall. If the wall is advancing into the fluid with a velocity U, then 
the mass flux across the interface (from Eq. 200) is pU. Apart from small compress- 
ibility effects, this is exactly the rate at which fluid mass is absorbed by the moving 
wall. For sliding motion, Eq. 200 correctly predicts no net mass transfer across the 
interface. 



4.2 

Particle motion 

An explicit update of the particle velocity 

U(t + h) = U(t) + — F(t) (202) 

m 

has been found to be unstable [99] unless the particle radius is large or the particle 
mass density is much higher than the surrounding fluid. In previous work [99] the 
instability was reduced, but not eliminated, by averaging the forces and torques over 
two successive time steps. Subsequently, an implicit update of the particle velocity 
was proposed [134] as a means of ensuring stability. A generalized version of that 
idea, which can be adapted to situations where two particles are in near contact, 
was developed in Ref. [104]. Here we sketch an elaboration of this idea, which is 
consistent with a Trotter decomposition of the Liouville operator [135-139]. We will 
only consider the update of the position and linear velocity explicitly; the extension 
to rotational motion is straightforward [104]. 

The equations of motion for the suspended particles are written as 

R, = Ui, (203) 
™U, = F? (Ri , Ui) + F c t (R N ), (204) 

where we have separated the forces into a hydrodynamic component F i , which de- 
pends on the particle position and velocity, and a conservative force F?, which de- 
pends on the positions of all particles. The hydrodynamic force depends on the fluid 
degrees of freedom as well, but these remain unchanged during the particle update 
and need not be considered as dynamical variables here. 

A second-order Trotter decomposition [135-139] breaks the update of a single 
time step into three independent components: a half time step update of the positions 
at constant velocity, a full time step update of the velocities with fixed positions, and 
a further half time step update of the positions using the new velocities: 

Ri(i + \h) = Ri(t) + |u<(t), (205) 



m 



F^(R 4 (t + \h), Ui) + F?(R N (t + \h)) 

h 
2 



(206) 



Ri(t + h)= Ri(t + \h) + -Ui{t + h). (207) 
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In the absence of velocity-dependent forces this is just the Verlet scheme, but the 
solid-fluid boundary conditions (Eq. 200) introduce a hydrodynamic force that de- 
pends linearly on the particle velocity [104, 134], 

FftRi, Ui) - F h (R 4 ) - C(Ri) • Ui. (208) 

The velocity independent force is calculated at the half-time step 

Fj(R i (t+i/ l )) = ^-^2<(r,t)c 6 , (209) 

b 

where the sum is over all the boundary nodes, b, describing the particle surface and Cb 
points towards the particle center. The location of the boundary nodes is determined 
by the particle coordinates Rj(t + |ft), which should be evaluated at the half time 
step as indicated. The post-collision populations, n£, are calculated at time t but 
arrive at the boundary nodes at the half time step also. The components of the matrix 

C(R l (t + ^h)) = ^J2 aCbC b C b ( 21 °) 

c s n b 

are high-frequency friction coefficients, which describe the instantaneous force on a 
particle in response to a sudden change in velocity. Complete expressions, including 
rotation, are given in Ref. [104]. 

The LB fluid and the solid particles are coupled by an instantaneous momentum 
transfer at the half time step, which is therefore presumed to be conservative: 

Ui(t + \h) = Ui(t) + ^jF?(R w (i + \h)), (211) 

U*(i + \h) = Ui(t + \h) + ^F?(Rj(t + |ft), Ui(t + \h)), (212) 

Ui(t + ft) = U?(* + ±ft) + ^Fj(R w (t + ±ft)). (213) 

However it is not entirely clear what velocity should be used in Eq. 212: among 
the possibilities discussed in Ref. [104] are an explicit update Uj(t + |ft) = 
Uj(f + | ft), an implicit update XJi(t + |ft) = U*(i + |ft), and a semi-implicit 
update Ui(i + ^ft) = [U t (t + \h) + XJ*(t + \h)]/2. It has been pointed out [140] 
that, even for a Langevin equation with constant friction, there are deviations in the 
temperature for finite values of h. However the semi-implicit scheme satisfies the 
FDT exactly for constant friction. Here we will consider a different model for the 
velocity, assuming that the hydrodynamic force is distributed over the time step. For 
simplicity we consider a single component of the velocity, 

mil = -(U + Ffr + F c , (214) 

where £, Fq and F c are constant in this context. The solution of Eq. 214 over a time 
interval h is 

jph _i_ pc 

U(t + h) = U(t) exp(-a) + -2-^ [1 - exp(-a)] , (215) 
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where a = (h/m is the dimensionless time step. Equation 215 is stable for all 
values of a, satisfies the FDT exactly, and, when there is no conservative force, 
encompasses previous algorithms as limiting cases. Both explicit [98] and implicit 
[104, 134] schemes are consistent with an expansion of Eq. 215 to linear order in a, 
while the semi-implicit method [104] can be derived from a second order expansion 
in a. The steady-state velocity U(t + h) = U (t) satisfies the force balance F h + 
F c = exactly. This new result may lead to more accurate integration of the particle 
positions and velocities in the large a limit. 

To complete the update, the velocity Uj(f + ^h) is needed to calculate the mo- 
mentum transfer to the fluid (Eq. 200). An explicit update [98] can be done in a single 
pass since TJ i(t + ^h) = Uj(i) is already known, but an implicit or semi-implicit 
update requires two passes through the boundary nodes. The first pass is used to cal- 
culate Fq so that Eq. 212 can be solved for Uj(t + \h) [104]. This velocity is used 
to update the population densities in a second sweep through the boundary nodes. 
In the present case we calculate \Ji(t + \h) by enforcing consistency between the 
sequential update Eqs. 211-213 and Eq. 215: 

aUi{t+\h) = Ui(t) [1 - cM-a)} + (K + F c ) - i- *^-") ^ .(216) 

This ensures overall momentum conservation as before. 

When there are short-range conservative forces between the particles, the LB 
time step is frequently too large for accurate integration of the interparticle forces. 
The LB time step can be divided into an integer number of substeps, but the question 
then arises as to how to best incorporate the hydrodynamic forces, since F h should, 
in principle, be calculated at t + \h. One possibility is to use the fact that £ varies 
slowly with particle position and accept the small error associated with using R(t) 
rather than R(t + \h). Or this solution could be used as a predictor step for calculat- 
ing R(t + i/i), which could then be followed by one or more corrector cycles with 
increasingly more accurate calculations of £(t + \h). The corrector cycles should 
not involve a significant overhead since the boundary nodes would be largely the 
same from one cycle to the next, and the time-consuming lookup of LB population 
densities could be avoided. 

Although the momentum exchange between fluid and solid occurs instanta- 
neously at the half time step, in calculating U j (t + \ h) we have made the assumption 
that the hydrodynamic force is distributed over the time step. We actually attempted 
to derive an update for the velocity assuming that the hydrodynamic force acts over 
a very small fraction of the time step, but this has not lead to a sensible result as yet. 
It is not entirely clear if the assumption that the hydrodynamic force acts over the 
whole time step is valid, and does not, for example, produce an artificial dissipation. 
To resolve this question will require a detailed analysis of the fully coupled system, 
along the lines given in Sec. 4.5 for the simpler case of frictional coupling. A similar 
analysis for solid-fluid boundary conditions is an open area for further research. 
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4.3 

Surfaces near contact 

When two particle surfaces come within one grid spacing, fluid nodes are excluded 
from regions between the solid surfaces, leading to a loss of mass conservation. This 
happens because boundary updates at each link cause mass transfer across the solid- 
fluid interface, which is necessary to accommodate the discrete motion of the particle 
surface (see Sec. 4.1). The total mass transfer in or out of an isolated particle is 



regardless of the particle's size or shape. 

Although the sum J2b aCbc b is zero f° r an Y closed surface [104], when two parti- 
cles are close to contact some of the boundary nodes are missing and the surfaces are 
no longer closed. In this case AM ^ and mass conservation is no longer ensured. 
Two particles that remain in close proximity never reach a steady state, no matter 
how slowly they move, since fluid is constantly being added or removed, depending 
on the particle positions and velocities. If the two particles move as a rigid body mass 
conservation is restored, but in general this is not the case. The accumulation or loss 
of mass occurs slowly, and in many dynamical simulations it fluctuates with changing 
particle configuration but shows no long-term drift. However, we typically enforce 
mass conservation, particle-by-particle, by redistributing the excess mass among the 
boundary nodes [104]. An alternative idea is to ensure that there is always at least 
one fluid node in the gap between the particle surfaces. In dense suspensions it would 
be quite inaccurate to insert an artificial excluded volume around the particles, but a 
more promising idea is to cut back the particle surfaces along planes perpendicular to 
the line of centers [141, 142]. It remains to be seen if the hydrodynamic interactions 
retain the level of accuracy shown in Fig. 2. 

When two particles are in near contact, the fluid flow in the gap cannot be re- 
solved. For particle sizes that are typically used in multiparticle simulations (a < 56), 
the lubrication breakdown in the calculation of the hydrodynamic interaction occurs 
at gaps of the order of 0.1a. However, in some flows, notably the shearing of a dense 
suspension, qualitatively important physics occurs at smaller separations, typically 
down to 0.01a. Here we outline a method to implement lubrication corrections into 
a lattice-Boltzmann simulation. 

For particles close to contact, the lubrication force, torque, and stresslet can be 
calculated from a sum of pairwise-additive contributions [29], and if we consider 
only singular terms, they can be calculated from the particle velocities alone [143]. 
In lattice-Boltzmann simulations [100, 144] the calculated forces follow the Stokes 
flow results down to a fixed separation, approximately equal to the grid spacing b, and 
remain roughly constant thereafter (see Fig. 2). The simplest lubrication correction 
is to take the difference between the lubrication force at a gap s and the force at some 
cut off distance s c ; i.e. 




(217) 




(218) 
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F' = s > s c , (219) 

where U12 = Ui — U2, s = |R, 12 1 — a\ — a 2 is the gap between the two surfaces, and 
the unit vector R12 = Ri2/|R-i2|- Numerical tests of this procedure for the older 10- 
moment LB model are reported in Ref. [104]. Results for the MRT model are shown 
in Fig. 3, using a cutoff distance s c = 1.1b for the parallel component and s c = 0.7b 
for the perpendicular component. Even this simple form for the correction gives an 
accurate description of the lubrication regime, with the largest deviations occurring 
near the patch points. A more accurate correction can be obtained by calibrating each 
distance separately as in Stokesian dynamics and related methods [21, 145]. 




Fig. 3. Hydrodynamic interactions including lubrication, with particles of radius a = 2.5b. 
The solid symbols are the LB friction coefficients, ^ pU and £P er J> 5 for the relative motion of 
two spheres along the line of centers (left) and perpendicular to the line of centers (right). 
Results are compared with essentially exact results from a multipole code [21] in the same 
geometry (solid lines). 



4.4 

Improvements to the bounce-back boundary condition 

The bounce-back boundary condition remains the most popular choice for simula- 
tions of suspensions, because of its robustness and simplicity. The results in Figs. 2 
and 3 show that accurate hydrodynamic interactions, within 1-2%, can be achieved 
with quite small particles, particularly when combined with the MRT model. The 
reason that bounce-back works so well, despite being only first-order accurate, is 
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that the errors in the momentum transfer tend to cancel when averaged over a ran- 
dom sampling of boundary node positions [1 14]. In fact bounce-back can sometimes 
be more accurate than interpolation, where the errors, though locally smaller, do not 
cancel. 
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Fig. 4. Variation in friction coefficients with grid location for particles of radius a = 2.5b. 
The solid symbols are the variance in the LB friction coefficients, Zi< p " and A( perp for the 
relative motion of two spheres along the line of centers (left) and perpendicular to the line of 
centers (right). Results are a single standard deviation in the friction coefficients, calculated 
from 100 independent positions with respect to the grid. 



The most important deficiency of the bounce-back algorithm is the dependence 
of the force on the position of the nodes with respect to the grid. The results in Figs. 
2 and 3 are averages over 100 independent configurations, in which the relative po- 
sitions of the particles are the same but the pair is displaced randomly with respect 
to the underlying lattice. However, the variance in the friction for randomly sampled 
grid locations is small, typically of the order of 1 %, as can be seen in Fig. 4. Nev- 
ertheless there is a much larger fluctuation in the force around the particle surface, 
which is particularly problematic if the particles are deformable [146, 147]. Thus 
while the bounce-back method is quite accurate on average, locally the errors can 
be large. A detailed analytical and numerical critique of the bounce-back algorithm 
can be found in Ref. [148], together with an analysis of several of the modifications 
mentioned below. The most practical higher-order boundary conditions are adapted 
from the link bounce-back algorithm outlined in Sec. 4.1. 

More sophisticated boundary conditions have been developed using finite-volume 
methods [149, 150] and interpolation [151-153]. A simple, physically motivated in- 
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terpolation scheme has been proposed [151, 154], which both improves the accuracy 
of the bounce-back rule and is unconditionally stable for all boundary positions; the 
scheme has both linear and quadratic versions. A more general framework for this 
class of interpolation schemes has been extensively analyzed in a comprehensive and 
seminal paper [1 13]; the Multi-Reflection Rule proposed in Ref. [1 13] is the most ac- 
curate boundary condition yet discovered for lattice-Boltzmann methods. However 
interpolation requires additional fluid nodes in the gap between adjacent particle sur- 
faces. The bounce-back rule requires only one grid point between the surfaces but 
linear interpolation requires at least two grid points, while quadratic interpolation 
and multi-reflection require three. Recently, it was proposed that only the equilib- 
rium distribution needs to be interpolated [114]. Although this is more complex to 
implement than linear interpolation, it has the advantage that the velocity distribution 
at the boundary surface may be used to provide an additional interpolation point. In 
this way the span of fluid nodes can be reduced to that of the bounce-back rule, while 
obtaining second-order accuracy in the flow field. In conjunction with an appropri- 
ate choice of collision operator [113], the location of the hydrodynamic boundary 
remains independent of fluid viscosity, unlike the linear and quadratic interpolations 
[151]. For viscous fluids, where 7 S > 0, the equilibrium interpolation rule is more 
accurate than either linear or quadratic interpolation [1 14]. 

4.5 

Force coupling 

The force-coupling algorithm [121, 122] starts from a system of mass points which 
are coupled dissipatively to the hydrodynamic continuum. The particles are specified 
by positions r i7 momenta p i7 masses rrii, and phenomenological friction coefficients 
Fi. They interact via a potential V (r N ), giving rise to conservative forces = 
—dV/dri. The fluid exerts a drag force on each particle based on the difference 
between the particle velocity and the fluid velocity u, = u(rj), 



Momentum conservation requires that an equal and opposite force be applied to the 
fluid. Both discrete and continuous degrees of freedom are subject to Langevin noise 
in order to balance the frictional and viscous losses, and thereby keep the temperature 
constant. The algorithm can be applied to any Navier-Stokes solver, not just to LB 
models. For this reason, we will discuss the coupling within a (continuum) Navier- 
Stokes framework, with a general equation of state p(p). We use the abbreviations 
Va/s-yS for the viscosity tensor (Eq. 46), and 

^aP = -P<W + PU a Uf3 (221) 

for the inviscid momentum flux or Euler stress (Eq. 79). Since the fluid equations 
are solved on a grid, whereas the particles move continuously, it will be necessary to 
interpolate the flow field from nearby lattice sites to the particle positions [122]. 




(220) 
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The addition of a point force into the continuum fluid equations introduces a sin- 
gularity into the flow field, which causes both mathematical and numerical difficul- 
ties. On the other hand, the flow field around a finite-sized particle can be generated 
by a distributed force located entirely inside the particle [155, 156]. This flow field 
is everywhere finite, and the force density appearing in Eq. 17 can be written as 

f(r) = -^Ffz\(r, ri ), (222) 

i 

where A(r) is a weight function with compact support and normalization 

d 3 rA(r,Ti) = 1. (223) 



Compact support limits the set of nodes r to those in the vicinity of and ensures 
that the interactions remain local. Away from solid boundaries, translational invari- 
ance requires that 

A(r,r i )=A(r-r i ). (224) 

The function A(r, r,) plays a dual role, both interpolating the fluid velocity field 
to the particle position, 

ufc) - J d 3 r/\(r,r,)u(r), (225) 

and then redistributing the reactive force to the fluid, according to Eq. 222. Within the 
context of polymer simulations, A has been regarded as an interpolating function for 
point forces, but it can equally well be regarded as a model for a specific distributed 
force, contained within an envelope described by A(r — rj). The flow fields from 
a point force and a distributed force are similar at large distances from the source, 
but the distributed source has the advantage that the near field also corresponds to 
a physical system, namely finite-size particles. We will adopt the distributed source 
interpretation both here and in Sec. 4.6. 

The Langevin equations of motion for the coupled fluid-particle system are: 

d 1 

-t: r * = — Pi, (226) 

at rrii 

j tPl = FC + Ft + F{, (227) 

dtp + d a j a = 0, (228) 
d t j a + dpi: = dprjapysdjUs + f£ + d (i a s a(} , (229) 

where the force density applied to the fluid includes both dissipative and random 
forces, 

f /l (r) = -^(Ff + Ff)z\(r,r,). (230) 

i 

The Langevin noises for the particles and fluid, F{ and c^, satisfy the usual moment 
conditions: 
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X>=0, (231) 
a£,; = (). (232, 



(*) (*')) = 2k B Tr i 6 ij 5 aP 6 (t - t') , (233) 

(a f aff (r, t) at (S (r', f )) = 2k B T Va0 , (S 6 (r - r') 5 (t - f ) . (234) 
By construction, this coupling is local, and conserves both the total mass 

M = m + ( d 3 rp (235) 

i 

and the total momentum 

P = E + / d3r P u - ( 236 ) 



Galilean invariance is ensured by using velocity differences in the coupling between 
particles and fluid (Eq. 220). A finer point is that the interpolation uses u (and not j), 
so that the velocity field enters strictly linearly. We will prove that the fluctuation- 
dissipation theorem (FDT) holds for this coupled system, proceeding in three steps 
that successively take more terms into account. 

Let us look first at the conservative system where the particles and fluid are com- 
pletely decoupled: 

^-r, = — Pi , (237) 

at m,{ 

= (238) 

d t p + d a3a = 0, (239) 

d t j a + d P 7T^ = 0. (240) 

The dynamics of the particles and the Euler fluid can be described within the frame- 
work of Hamiltonian mechanics [157]. The Hamiltonians for the particles 

Pi 



and fluid, 

H f = J rf 3 rQpu 2 +e(^ , (242) 

are conserved quantities, with e(p) the internal energy density of the fluid. 

As a second step, we consider a system where particles and fluid are still decou- 
pled, but are subject to dissipation and noise: 

d 1 

—ri = — Pi , (243) 

at nii 
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^= F -S pi+F - (244) 

dtp + d a j a = 0, (245) 

dtja + dpiTap = dpriafj-fsd-fUs + dpa^p. (246) 

These Langevin equations are known to satisfy the FDT [15, 42, 158-160]. We 
briefly sketch the formalism used for the proof, since this will be needed for the 
final step in which we consider the fully coupled system. 

Instead of describing the stochastic dynamics via a Langevin equation, we use the 
Fokker-Planck equation, which is the evolution equation for the probability density 
in phase space. For an N -particle system, 

8 t P (v N , p N ) = (d + £ 2 + £ 3 ) P (v N , p N ) , (247) 

where r N , p N denote the positions and momenta of all N particles. The three oper- 
ators C\, £2, and £ 3 describe the Hamiltonian, frictional, and stochastic part of the 
dynamics; they can be found via the Kramers-Moyal expansion [15, 159]: 

A = -W^ + ^- F A (248) 



m l dpi 

d 

rrii dp t 



— ' m ■ rir\: 



Cs = k B T^2r~. (250) 

i * 

The FDT holds if the Boltzmann factor, exp (— Hp/ksT), is a stationary solution of 
the Fokker-Planck equation. Using (i = (fc^T) -1 to define the inverse temperature, 
we have 

£1 exp {-pHp) = (251) 

as a direct consequence of energy conservation in Hamiltonian systems. Further- 
more, the relation 

(£ 2 + £ 3 ) exp (-/?W P ) = (252) 

can be shown by direct differentiation. 

For the fluid system, the phase space comprises all possible configurations of the 
fields p(r), j(r), which we denote as [p] , [j]. The Fokker-Planck equation for the 
fluid degrees of freedom can be written as 

d t P([p] , [j]) = (£4 + £5 + £ 6 ) P ([p] , , (253) 
where £4, £5, and C§ describe the Hamiltonian, viscous, and stochastic components, 

£4 = j d 3 r (Jjdaja + -Jj-dp^ , (254) 
(pr—dpd-fUs, (255) 



ksTvap-yS J d 3 r J d 3 r' 



Or 13 dr' s 



(256) 
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and 5 . . . /S . . . represents a functional derivative [161]. Replacing d/dr' s in the last 
equation with —d/drg enables integration over r': 

£e = -ksTriapys J d 3 r-^-df!d 7 -^-, (257) 

where we have exploited the symmetry of the viscosity tensor with respect to the 
indexes 7 and 5. Functional differentiation of the Boltzmann factor with respect to j, 

— exp(-pH f ) = -l3usexp(-pHf), (258) 
ojs 

then shows that 

(£ 5 + £ 6 )cxp(-/3W / ) =0. (259) 
Finally, the relation 

£ 4 exp(-/?W / ) = (260) 

follows from energy conservation in Hamiltonian dynamics. 

We now turn to the coupled system, with Hamiltonian H = H p + Hf. The 
Fokker-Planck equation in the full phase space reads 

d t P (r N , p N , [p] , [j]) = (j2 P P^' W ' K) ' < 261 > 

with the operators £7 — £ w to describe the coupling in the equations of motion: 

d 

'' dp ia 



£7 = - ^2 r iJT— U iai ( 262 ) 



£ 8 = - V P j d 3 rA (r, Ti ) — ^ (—Pia - Via) , (263) 
£, = k B TY j r i JfiTA^^jMAV^)^ (264) 

£ W = -^BT^r.jL /^(r.rO (265) 

The coupling of the fluid velocity to the particles is described by £7, while £ 8 de- 
scribes the drag on the fluid. The stochastic contributions include fluid-fluid correla- 
tions via £9, and fluid-particle cross correlations via £10. 

It should be noted that u^, being the result of the interpolation, depends on the 
fields [p] and [j], so that 8ui/Sj a (r) is nonzero. Hence, the corresponding operators 
in £ 8 do not commute. Explicit functional differentiation shows that Eq. 258 holds 
in an analogous way for the interpolated velocity, 

J d 3 vA (r, r< ) j^j- exp (-0H) = -pu ia exp (-0H) . (266) 

Explicit calculations, as outlined for the uncoupled system, show that 
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(d + d + d + C w ) cxp (-/3H) = 0, (267) 
which implies the exact FDT for the fully coupled system, 

(^2c)jcxp(-f3H) = 0. (268) 

An important consequence of this result is that a consistent simulation needs to ther- 
malize both fluid and particle degrees of freedom; any other choice will violate the 
FDT. 

The Langevin integrator for the particles is constructed in much the same way as 
the velocity Verlet algorithm for Molecular Dynamics. Although a Langevin analog 
to the Verlet algorithm has been known for some time [162], straightforward deriva- 
tions have become available only recently, by applying operator-splitting techniques 
that were previously limited to Hamiltonian systems [163]. We employ a second- 
order integrator [135-137], which reduces to the velocity Verlet scheme in the limit 
of vanishing friction. Higher-order schemes are known [164], but they are consider- 
ably more complicated. Specifically, we approximately integrate the equations 

4 r < = —Pi. ( 269 > 

at rrii 

l Pi = F? + F? + F{, (270) 

assuming that the fluid velocity is constant over a time step h. This corresponds 
to the Fokker-Planck equation for the particles 

d t P (r^p", t) = (d + d) P (r W , P V) (271) 
with 

4- = -E/r---. (272) 

= - E w, • F < + • " u ') + k ^ r 'M ,273) 

The formal solution 

P (r N ,p N ,h)= cxp [(d + d) h ] p - P N > 0) (274) 
is approximated by a second-order Trotter decomposition, 

cxp [(d + d) h \ = CX P {dh/2) cxp (C p h) cxp (C r h/2) + 0(h 3 ). (275) 

The operator splitting implies the following algorithm: a half time step update of the 
coordinates with constant momenta, 

r i (t + h/2) = v i (t) + ^^-, (276) 

followed by a full time step momentum update, with constant coordinates, 

Pi (t + h) = cVWpi® + C\ 2 \h) [F 4 C + r^] + C\ 3 \h)0i, (277) 
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and finally another half time step coordinate update, with constant momenta, 

v % (t + h) = v l {t + h/2)+ hlpl{t + h) 
The coefficients in Eq. 277 are 



TO, 



Cf ) (/ i )=exp 



Cf\h) = 



Cf\h) 



1 — exp I ——h 

' mi 



imik B T 



1 — cxp 



mi 



(278) 

(279) 
(280) 

(281) 



where 0i a are Gaussian random variables with zero mean and unit variance. Equa- 
tion 277 is the exact solution of the momentum update, since Eq. 270 is a linear 
Langevin equation describing Brownian motion in a harmonic potential [158]. Thus 
the only source of error in integrating the particle motion is derived from the Trotter 
decomposition itself, Eq. 275. Nevertheless, it is important to limit the range of ran- 
dom numbers, to ensure that very large steps do not occasionally occur. It is therefore 
both desirable and more efficient to use distributions of random variates with finite 
range, which reproduce the Gaussian moments up to a certain order; in the present 
case, 0{h 2 ) accuracy requires that moments up to the fourth cumulant are correct. 
One possible choice is 



p(0) = \s (0) + ±s (o - Vs) + (e + Vz) 



(282) 



After the update of the particle momenta, the fluid force density, f h (r) (Eq. 230), 
is distributed to the surrounding lattice sites. After one or more particle updates, the 
fluid variables are updated for a single LB time step, with the external forces being 
taken into account in the collision operator. This scheme is probably only first order 
accurate overall. It remains a challenge for the future to develop a unified frame- 
work to describe the fully coupled system and analyze its convergence properties; 
the algorithm could then perhaps be improved in a systematic fashion. 

The input friction coefficient is not the same as the long-time friction coefficient, 
which is measured by the ratio of the particle velocity to the applied force [122]. 
Consider an isolated particle with "bare" (or input) friction coefficient r dragged 
through the fluid by a constant force F, resulting in a steady particle velocity U. The 
force balance requires that 

F = -F d = r(U - uo), (283) 

where u is the fluid velocity at the particle center, r , 



u = 



J d 3 rZ\(r,r )u(r). 



(284) 



In the absence of thermal fluctuations, the deterministic fluid velocity field can be 
calculated in the Stokes flow approximation using the Green's function appropriate 
to the boundary conditions [22], 
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u(r) = J d 3 r'T(r,r') • FA(t',t q ). (285) 



In an unbounded fluid, the Green's function reduces to the Oseen tensor [22], 
T(r,r') = O(r-r'), with 

1 rr\ 

&TTTir V r 2 

but Green's functions are also known for periodic boundary conditions [165] and 
planar boundaries [38, 39] as well. Combining Eqs. 283, 284 and 285, 

uo = T av ■ F = n av • (U - u ) , (287) 

where T av is the Green's function averaged over the particle envelope, 

T^r-o) = Jd 3 r J d 3 r'Z\(r,r )Zi(r',r )T(r,r'). (288) 

In a system with translational invariance, T av is independent of r and propor- 
tional to the unit tensor by symmetry. Then from Eq. 287 

F = — - r U, (289) 

1 + flool 

where ^ — T™/3 accounts for the renormalization of r. If the size of the particle 
a is associated with the range of interaction of A, then dimensional analysis of Eq. 
288 suggests that fi^ <~ [r]a)~ 1 . The effective friction coefficient, defined via F = 
r e f fXJ, is therefore diminished by the flow field induced by the applied force, 

! ^+Moo- (290) 



r eff r 

This relation has been verified numerically by extensive computer experiments 
[122, 140]. In the strong coupling limit, r — > oo, the effective friction saturates 
to the limiting value r e ff = [i^. This suggests assigning an effective radius to the 
interpolating function, 

- = 6tt?7 Moo . (291) 
a 

However, for smaller values of the input friction, the effective particle size is given 
by 

I = 67n? ( I + Moo ) = 1 + J_ ; (292) 
a \T ) a gb 

where a = Fj 6iri] is the input particle radius and gb = (6tti]^i 00 )^ 1 depends on the 
interpolating function. The interesting physical parameter is r e ff, which describes 
the long-time behavior of the coupled system. Thus, to approach the continuum 
limit, one should keep r e f / constant as the lattice spacing is decreased, and change 
the bare coupling r as necessary. 

It is not yet known how T av (ro) behaves in the vicinity of a solid boundary, 
when translational invariance is broken. The compact support of the weighting func- 
tion A suggests that T av is a local correction and therefore largely independent of 
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macroscopic boundary conditions. Numerical simulations with periodic boundary 
conditions (Sec. 4.6) show that g is independent of system size and fluid viscosity. 
The weak system-size dependence reported in Ref. [140] is entirely accounted for 
by the difference between the periodic Green's function [165] and the Oseen tensor. 
Thus in a periodic unit cell of length L, Eq. 292 requires a correction of order 1 /L 
[19], 

1 1 2.84 1 



gb a L ao 



(293) 



4.6 

Interpolating functions 

In this section we consider translationally invariant interpolating functions, A(r, Tq) 
= A(r — r ) in more detail. For a discrete lattice, the interpolation procedure (c.f. 
Eq. 225) reads 

u(r-o) =^Z\(r-r )u(r), (294) 

r 

where ro is the position of the particle, and r denotes the lattice sites. The normal- 
ization condition 

J2 ^(r - r ) = 1 (295) 

r 

must hold for all particle positions, r , in order for the conservation laws to be sat- 
isfied exactly. We will assume the analysis of Sec. 4.5 carries over to the discrete 
system with no more than second-order discretization errors. Proof of this assump- 
tion remains for future work; here we numerically compare various choices of A. 

Previous work, incorporating force coupling into spectral codes, used an isotropic 
Gaussian distribution for A(r) [155], but this is not commensurate with cubic lattice 
symmetry. Thus we take A as a product of one-dimensional functions [129] 

A(x,y,z) = 4>(^)4>(l)<t>(l)- (296) 
We first consider the two-point linear interpolating polynomial 

fc(«) = | | U |> 1; (297) 

which satisfies the following moment conditions for all real-valued u and integer j: 

X>(«-J) = 1 . ( 298 > 

j 

J2j<f>(u-j)=u. (299) 

3 

It exactly conserves momentum and angular momentum of the particle and fluid. 
However, <p 2 violates the condition 



54 



Burkhard Diinweg and Anthony J. C. Ladd 



»(«) 



J2<P 2 (u-j)=C, (300) 

3 

where C is a constant, independent of u. The importance of this condition is ex- 
plained in Ref. [129]. 

The conditions in Eqs. 298-300 can be satisfied by a three-point interpolation 
function, 

i(l + Vl-3u 2 ) o<M<± 
\ [5 - 3\u\ - y/-2 + 6\u\ - 3u 2 ) \ < \u\ < § (301) 

| < |u|. 

An important property that emerges from fa is that the first derivative, <f>' 3 (u), is 
continuous throughout the whole domain of u. This ensures that the velocity field 
varies smoothly across the grid, with a continuous spatial derivative Vu. By contrast, 
linear interpolation leads to a continuous velocity but discontinuous derivatives. 

In order to test the various interpolation schemes we have determined the settling 
velocity of a single particle in a periodic unit cell. A small force was applied to 
the particle and a compensating pressure gradient (or uniform force density) was 
added to the fluid, so that the net force on the system was zero. The steady-state 
particle velocity was determined, without allowing the particle to move on the grid 
[99]. This procedure is valid in Stokes flow, where an arbitrarily small velocity may 
be assumed, and gives a clean result for the variation in settling velocity with grid 
position. We used Eq. 293 to convert the measured mobility to a single parameter 
g, which does not depend on system size (L) or fluid viscosity (77). The mobility of 
the particle is reduced by the periodic images [165], but the correction in Eq. 293 
accounts for the effects of the periodic boundaries quantitatively [99, 120, 166]. In 
simulations of polymer solutions an average of g over all grid positions is used. 

The smoother velocity field derived from three-point interpolation means that the 
particle velocity is less dependent on the underlying grid than with linear interpola- 
tion, as can be seen in Fig. 5. Here we show the variation in effective particle size, 
as determined by the parameter g (Eq. 293), for linear (two point) interpolation (Eq. 
297), three-point interpolation (Eq. 301), and the four-point interpolation, 



f(3-2M + v /l + 4H-4u 2 ) 0<|it|<l 
Mu) = I 1(5- 2\u\ - V-7+ 12M -4m 2 ) 1 < |u| < 2 

2<|u|, 



(302) 



that is commonly used in immersed boundary methods [129]. As was recently no- 
ticed [130], the four-point interpolation leads to a much smaller variation in effective 
friction than linear interpolation. The parameter g varies with grid position by up to 
20% in the case of linear interpolation, but by less than 1 % with four-point interpo- 
lation. On the other hand, linear interpolation requires an envelope volume of 8 grid 
points, while the four-point scheme requires 64 grid points. Away from the strong- 
coupling limit, the grid dependence of the settling velocity is reduced since there is 
a non-negligible input mobility apart from the lattice contribution. 
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Fig. 5. Variation in settling velocity with grid location. The effective hydrodynamic radius 
was determined from U = F/Qnrja and converted to g using Eq. 293; g was found to be 
independent of r] and L, as expected. Results are shown at 56 different grid positions (labeled 
by the index n), systematically varying the coordinates in steps of 0.16. Particles with input 
radius do = b were placed at coordinates (ib/10, jb/10, kb/10), with 0<i<j<fc<5. 
Results are shown for the two-point (circles), three-point (squares) and four-point (diamonds) 
interpolation schemes. 

Four-point interpolation is only necessary when using centered-difference ap- 
proximations to the velocity and pressure fields [129], a situation that does not arise 
in LB simulations. We see that the three-point scheme is also much smoother than 
linear interpolation, with about a 3% variation in g. It is not as smooth as the four- 
point interpolation, but requires only 27 grid points. Furthermore, the smaller span of 
nodes means that the boundary surface is more tightly localized, and in fact the hy- 
drodynamic interactions obtained with three-point interpolations are just as accurate 
as those obtained with four-point interpolation, as shown below. 

An important test of the force coupling scheme is its ability to represent the 
hydrodynamic interactions between two spherical particles. As an example of the 
accuracy of the different interpolation schemes, in Fig. 6 we show the hydrody- 
namic interactions between two spheres moving along the line of centers. A small 
force is applied to sphere 1, in the direction of the vector between 1 and 2, and 
the velocity of sphere 2 is determined. From this we can calculate the hydrody- 
namic mobility fjffi ■ The results are normalized by the mobility of the isolated sphere 
/io = (6777701) _1 . Results were obtained for the two-point, three-point, and 4-point 
interpolation schemes, using a source particle placed on a grid point (0, 0, 0) in one 
instance and in the center of the voxel (6/2, 6/2, 6/2) in the other. The simulations 
were carried out in a periodic unit cell, with 20 grid points in each direction. Results 
are compared with a spectral solution of the Stokes equations in a periodic geom- 
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etry, assuming the force density on the particle surface is constant. For an isolated 
pair of particles this level of approximation includes the Oseen interaction and the 
Faxen correction; it corresponds to the Rotne-Prager (RP) interaction [23] used in 
most Brownian dynamics simulations of hydrodynamically interacting particles: 

„2 



U 2 



T(Rl2) + i2^^ ( ^ 2l ~ 3Rl2Rl2) 



•Fi. (303) 



The periodic RP tensor can be calculated by Ewald summation [167] or by direct 
summation of Fourier components, which converges if used in conjunction with 
finite-volume sources [17, 168]. 

The results for the two-point interpolation show a significant dependence on the 
exact grid position when the particles are close to each other, r < 4a, but the results 
for the higher-order schemes are essentially independent of the grid. The three-point 
integration scheme is of comparable accuracy to the four-point scheme but requires 
less than half the number of grid points. It would seem to be the best choice for ap- 
plications even though the particle motion is not quite as smooth. When the particles 
are widely separated, the simulations match almost perfectly with both the Oseen 
and RP solutions for the same periodic geometry; the typical errors are of the order 
of 0.1% of the Stokes velocity. 

When the spheres are closer together, r < 3a, then the simulated hydrodynamic 
interactions match the Rotne-Prager interaction rather better than the Oseen interac- 
tion. This confirms that the weight function does make the particles behave as volume 
sources, rather than points. The best fit between simulation results and Stokes flow is 
obtained for an effective particle radius that is roughly 0.33w, where w is the range 
of the weight function. So for two-point interpolation (w = 26) the effective size is 
about 0.76, for three point interpolation (w = 36) it is about 1.06 and for four-point 
interpolation (w = 46) it is about 1.36. The actual values used in Fig. 6 are 0.86, 
1.06 and 1.26 respectively. The optimal hydrodynamic radius is quite large, a ~ 6, 
corresponding to strong coupling, ao ~ 56. It remains an open question whether it is 
practical or desirable to run the fluctuating simulations with such large input friction. 



5 

Applications with hydrodynamic interactions 

In this section we will discuss applications of the LB method to simulations of soft 
matter. We will briefly summarize some of our published work in this area, with the 
aim of indicating the breadth of possible applications of the method. Results will be 
summarized from simulations that cover a wide range of experimental length and 
time scales-from nanometers to millimeters and from nanoseconds to seconds. 



5.1 

Short-time diffusion of colloids 



The first application of the fluctuating LB model was to the short-time diffusion of 
hard-sphere colloids. At the time, a new experimental technique-Diffusing Wave 
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Fig. 6. Hydrodynamic interactions between a pair of spherical particles using the force- 
coupling method. The normalized mobility /J^/^io is plotted at various separations r. We 
fit the effective hydrodynamic radius of each interpolating function to numerical solutions of 
the Stokes equations with a uniform force density on the sphere. These results correspond to 
a Rotne-Prager description of the hydrodynamic mobility and do not include lubrication. We 
show results at two different grid locations, (0, 0, 0) (circles) and (6/2, 6/2, 6/2) (squares), 
for two-point (left), three-point (center) and four-point (right) interpolation schemes. 

Spectroscopy (DWS) [169]-enabled the study of the dynamics of colloids on time 
scales of a few nanoseconds. In general, the diffusion of a colloidal particle is a 
Markov process, but at such short times the developing hydrodynamic flow field 
gives rise to additional long-range correlations, analogous to the "long-time tails" 
in molecular dynamics [9]. Although the existence of long-time tails had been es- 
tablished theoretically [170, 171] and by molecular dynamics simulations [9], these 
experiments marked the first direct observation of correlated hydrodynamic fluctua- 
tions. Brownian and Stokesian dynamics both neglect long-range dynamic correla- 
tions, using instead the Stokes-flow approximation, which is typically only valid on 
time scales longer then one microsecond. 

Long-time tails occur naturally in the dynamics of lattice-gas models of col- 
loidal suspensions [173] and even of the lattice gases themselves [174]. It might be 
supposed that such correlations would be absent in a Boltzmann-level model, due 
to the Stosszahlansatz closure assumption. The fluctuating LB model described in 
Sec. 3 does not have any long-time tails in the stress autocorrelation functions, but 
mode coupling between the diffusion of fluid momentum and the diffusion of the 
colloidal particle does lead to an algebraic decay of the velocity correlation func- 
tion of a suspended sphere [175]. In these simulations the particle-fluid coupling 
was implemented via the link-bounce-back (BB) algorithm [98, 120] described in 
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Fig. 7. Decay of translational (U) and rotational (J?) velocity correlations of a suspended 
sphere. The time-dependent velocities of the sphere are shown as solid symbols; the relaxation 
of the corresponding velocity autocorrelation functions are shown as open symbols (with sta- 
tistical error bars). A sufficiently large fluid volume was used so that the periodic boundary 
conditions had no effect on the numerical results for times up to t = 1000 in lattice units 
(h — b — 1). The solid lines are theoretical results, obtained by an inverse Laplace transform 
of the frequency dependent friction coefficients [172] of a sphere of appropriate size (a = 2.6) 
and mass (p s /p = 12); the kinematic viscosity of the pure fluid rj kin = 1/6. 

Sec. 4.1. Figure 7 shows the decay of translational and rotational velocity from two 
different types of computer experiment. In one case an initial velocity is imposed, 
which decays away due to viscous dissipation, and in the other the particle is set 
in motion by stress fluctuations in the fluid. Figure 7 shows that, within statistical 
errors, the normalized velocity correlation functions are identical to the steady de- 
cay of the translational and rotational velocities of the sphere; thus our simulations 
satisfy the fluctuation-dissipation theorem. Moreover the simulations agree almost 
perfectly with theoretical results derived from the frequency-dependent friction co- 
efficients [172], even though there are no adjustable parameters in these compar- 
isons; thus we see that the fluctuating lattice-Boltzmann equation can account for the 
hydrodynamic memory effects that lead to long-time tails [9]. 

The simulations were also used to measure self diffusion in dense colloidal sus- 
pensions, up to a solids volume fraction of 45%. The simulation data, shown in 
Fig. 8, exhibits the same scaling with amplitude and time found in the DWS ex- 
periments [176]. In Fig. 8 the amplitude of the mean-square displacement has been 
normalized by its limiting value 6D s t, where D s is the short-time self-diffusion co- 
efficient. D s {(j)) is a monotonically decreasing function of concentration, because 
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Fig. 8. Scaled mean-square displacement (AR 2 (t)^ /6D s t at short times, vs. reduced time 
t/r. Simulation results for 128 spheres (solid symbols) are shown at packing fractions (f> of 
5%, 25% and 45%; the solid line is the isolated sphere result. The suspension viscosity at these 
packing fractions is 1.14r?o, 2.17?7o, and 5. 6770 respectively, where r)o is the viscosity of the 
pure fluid. 



neighboring particles increasingly restrict the hydrodynamic flow field generated by 
the diffusing particle. Although the colloidal particles are freely moving in the fluid, 
the no-slip boundary condition induces stresslets and higher force multipoles on the 
particle surfaces; D s (4>) is one average measure of these hydrodynamic interactions. 
The normalized mean-square displacement has a single relaxation time, so that when 
the time axis is scaled by r all the data collapse onto a single curve, which is the same 
as for an isolated sphere. The relaxation time r is the viscous diffusion time pa 2 /rj, 
of a single particle in a fluid of viscosity 77, where i]((f)) is numerically similar to the 
high-frequency viscosity of the suspension [19]. This observation, which is in line 
with experimental measurements [176, 177], suggests that the short- time diffusion is 
essentially mean field like. 

5.2 

Dynamic scaling in polymer solutions 

The classical theory [6] of the equilibrium dynamics of polymer chains in solution 
is the Zimm model [178], which considers a single flexible chain in a good solvent, 
such that its conformations are given by a random coil with excluded volume seg- 
ments: 
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R~bN». (304) 

Here, R is the size of the coil, measured in terms of the gyration radius or the end-to- 
end distance, while N denotes the number of monomers in the chain or the degree of 
polymerization, and b is the monomer size. Long-range interactions like electrostat- 
ics, or effects of poor solvent quality, are not considered. Furthermore, the solution 
is considered to be dilute, such that the chains do not overlap and a single-chain 
picture is sufficient. In other words, the standard Zimm model applies in the upper 
left corner of the generic phase diagram given in Fig. 9. Here the exponent v takes 
the value v w 0.59 in three dimensions, because the excluded volume interaction 
leads to swelling of the chain when compared to an ideal random coil (y = 1/2). A 
polymer with excluded-volume is thus a self-similar random fractal with dimension 
l/v. 

The Zimm model is based on the Rouse model [6, 179], but includes long- 
range hydrodynamic interactions between the segments. Both models predict self- 
similarity, not only with respect to space, but also with respect to time. Therefore the 
dynamics is conveniently described in terms of an exponent z, connecting the chain 
relaxation time tr with the size of the coil R: 

t R (xR z . (305) 

The internal degrees of freedom completely re-organize on a time scale tr, leading 
to statistically independent conformations. This is also the time that the chain needs 
to diffuse through a distance equal to its own size: 

D cm TR ~ R 2 , (306) 

where D cm is the translational diffusion coefficient, describing the center-of-mass 
motion. These two relations can be combined to determine the scaling of D cm with 
chain size: 

D cm oc R 2 - z . (307) 

In the Zimm model, the hydrodynamic interactions result in strongly correlated mo- 
tions, such that the coil as a whole behaves like a Stokes sphere, D cm oc i2 _1 , or 

z = 3. (308) 

The Rouse model neglects hydrodynamic interactions and the monomer friction co- 
efficients add up to give a total friction coefficient that is linearly proportional to 7Y. 
Since D cm oc iV" 1 oc Rr x l v , 

z = 2+-, (309) 
v 

corresponding to slower dynamics. 

Self-similarity implies that the relaxation of the internal degrees of freedom also 
scales with the exponent z on time scales n, <C t <C tr, where t& is the relaxation 
time on the monomer scale b. In the space-time window b -C I -C R, n -C t <C tr, 
there is a scaling of the mean square displacement of a monomer, 

(Ar 2 ) oc t 2/z , (310) 
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Fig. 9. Phase diagram of a polymer solution, in the c — T plane, where c is the monomer 
concentration and T is the temperature (parameterizing solvent quality). The static properties 
are characterized by scaling laws which describe the dependence of the chain size R (gyra- 
tion radius or end-to-end-distance for example) on the degree of polymerization N. In the 
dilute limit (c — > 0) the so-called theta transition occurs, where at T — 0, single isolated 
chains collapse from a swollen random coil to a compact globule. For finite chain length N, 
this transition is "smeared out" over a temperature region AT oc N^ 1 ^ 2 , in which the chain 
conformations are Gaussian. Below O, there is phase coexistence between a "gas" of globules 
and a "liquid" of strongly interpenetrating Gaussian chains. The corresponding critical point 
occurs at a very low concentration, c c oc N^ 1 ^ 2 , and in the vicinity of 0, — T c oc N^ 1 ^ 2 . 
The crossover region, which connects the regime of swollen isolated coils with that of the 
concentrated (Gaussian) solution at high temperatures, is called the semidilute regime. The 
dynamics is characterized by the Zimm model in the dilute limit where hydrodynamic inter- 
actions are important, and by the Rouse model for dense systems where they are screened. For 
very dense systems or sufficiently long chains, where curvilinear motion dominates, the Rouse 
model must be replaced by the reptation model (or the crossover behavior between these two 
cases). The Rouse and Zimm models are briefly described in the text. 
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while the dynamic structure factor of a single chain of TV monomers 



^,*) = ^(Eexp(ik.(r i (i)-r j (0))) 




ij = l 



(311) 



scales as 



S(M) = fc- 1 /"/ (fc 2 * 2A ) • 



(312) 



The Zimm model applies to dilute solutions, and, therefore, to the dynamics of 
a single solvated chain. It has become a benchmark system, used to test the validity 
of mesoscopic simulation methods. A single chain, modeled by bead-spring interac- 
tions, coupled to a surrounding solvent to account for hydrodynamic interactions, has 
been successfully simulated via (i) Molecular Dynamics [180-182], (ii) Dissipative 
Particle Dynamics [183, 184], Multi-Particle Collision Dynamics [185, 186], and 
by Lattice Boltzmann [122], applying the dissipative coupling described in Sec. 4.5. 
These studies are nowadays all sufficiently accurate to be able to clearly distinguish 
between Rouse and Zimm scaling. However, a more demanding goal is to verify not 
only the exponent, but also the prefactor of the dynamic scaling law. 

The Kirkwood approximation to the diffusion constant, 



(r is the monomer friction coefficient) can be calculated from a conformational av- 
erage of the hydrodynamic radius R H , 



Highly accurate results for a single chain in a structureless solvent have been ob- 
tained by Monte Carlo methods [187, 188]. However, a naive comparison will fail 
badly. The expression for Rh (Eq. 314) assumes an infinite system, but in a simula- 
tion the system is confined to a periodic unit cell, which is typically not substantially 
larger than the size of the coil. The effects of periodic boundaries can be accounted 
for quantitatively, by replacing the Oseen tensor with an Ewald sum that includes 
the hydrodynamic interactions with the periodic images [167]. The consequences of 
this have been worked out in detail for polymers [182] and colloids [19]. The main 
result is that the hydrodynamic radius must be replaced by a system-size dependent 
effective hydrodynamic radius; the leading-order correction is proportional to R/L, 
where L is the linear dimension of the periodic simulation cell. Interestingly, inter- 
nal modes, such as Rouse modes [6], where the motion of the center of mass has 
been subtracted, have a much weaker finite-size effect, which scales as L~ 3 , cor- 
responding to a dipolar hydrodynamic interaction with the periodic images [122]. 
Taking the finite-size effects into account, the predictions of the Zimm model are 
nicely confirmed. 

However, the Zimm model is no longer valid as soon as the chains start to over- 
lap. Here a double screening mechanism sets in: (i) Screening of excluded volume 
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Fig. 10. Scaling of the single-chain dynamic structure factor data, showing both Rouse and 
Zimm scaling [43]. The wave number k has been restricted such that only length scales above 
the blob size are probed < 1), while the size of the polymer chain as a whole does not yet 
matter (kRc > 1). The data are labeled according to the time regimes; solid symbols refer to 
the short-time regime below the blob relaxation time, t < T£, while open symbols are for later 
times t > T£. The upper curve is for Rouse scaling (z = 4) and the lower curve for Zimm 
scaling (z = 3). One sees that Zimm scaling works better in the short-time regime, while 
Rouse scaling holds for later times. 



interactions (Flory screening). In a dense melt, the chain conformations are not those 
of a self-avoiding walk, but rather those of a random walk (u = 1/2). Essentially, 
this is an entropic packing effect: A swollen coil would take too much configura- 
tion space from the surrounding chains. This effect can be understood in terms of 
a self-consistent mean-field theory, which is expected to work well for dense sys- 
tems where density fluctuations are suppressed [7]. (ii) Screening of hydrodynamic 
interactions. In dense melts, the dynamics is not Zimm-like, but rather Rouse-like, 
or governed by reptation [6]. Reptation occurs for long chains in dense systems, 
where topological constraints enforce an essentially curvilinear motion. We will not 
be concerned with these latter effects, but rather with the mechanism which leads 
to the suppression of hydrodynamic interactions. Based upon the results of com- 
puter simulations [43], we were able to develop a simple picture, which essentially 
confirmed the previous work by de Gennes [189], and completed it. The basic mech- 
anism is chain-chain collisions. A monomer encountering another chain will deform 
it elastically, inducing a stress along the polymer backbone instead of propagating 
the signal into the surroundings. Since the chain arrangements are random, the fluid 
momentum is also randomized, such that momentum correlations (or hydrodynamic 
interactions) are destroyed. 
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The crossover region between dilute and dense systems is called the semidilute 
regime. A semidilute solution is characterized by strongly overlapping chains which 
are however so long and so dilute that the monomer concentration can still be con- 
sidered as vanishingly small. Apart from b and R, there is now a third important 
length scale, the "blob size" £ with b <C £ <C R. Essentially, £ is the length scale on 
which interactions with the surrounding chains become important; this length scale 
controls the crossover from dilute to dense behavior. The chain conformations are 
characterized by v = 0.59 on length scales much smaller than £, while on length 
scales substantially above £ the exponent is v = 1/2. The challenge for computer 
simulations is that both behaviors need to be resolved simultaneously, which is only 
possible for N > 10 3 . Roughly 30 monomers are needed to resolve the random frac- 
tal structure within the blob, while another 30 blobs per chain are needed to observe 
the random walk regime. Furthermore, a many-chain system should be run without 
self-overlaps, and this leads to the conclusion [43] that the smallest system to simu- 
late semidilute dynamics contains roughly 5 x 10 4 monomers and 5 x 10 5 LB lattice 
sites. 

The picture which emerges from these simulations [43] can be summarized as 
follows. Initially, the dynamics is Zimm-like, even for length scales beyond the blob 
size. The reason is that hydrodynamic signals can spread easily throughout the sys- 
tem, and just drag the chains with them. This continues until chain-chain collisions 
start to play a role. The relevant time scale is the blob relaxation time oc £ 3 , i.e. 
the time a blob needs to move its own size. From then on, the screening mecha- 
nism described above becomes important, and the dynamics is Rouse-like. This is 
only observable on length scales beyond the blob size, since on smaller scales all 
dynamic correlations have decayed already. This is nicely borne out by single-chain 
dynamic structure factor data (see Fig. 10), and explains the previous observation of 
"incomplete screening" [190] in a straightforward and natural way. 

5.3 

Polymer migration in confined geometries 

Flexible polymers in a pressure-driven flow field migrate towards the center of the 
channel, because of hydrodynamic interactions. The local shear rate stretches the 
polymer and the resulting tension in the chain generates an additional flow field 
around the polymer. This flow field becomes asymmetric near a no-slip boundary and 
results in a net drift towards the center of the channel [191-193]. Recent simulations 
[193, 194] show that hydrodynamic lift is the dominant migration mechanism in 
pressure-driven flow, rather than spatial gradients in shear rate. Recently, we used 
numerical simulations to investigate a flexible polymer driven by a combination of 
fluid flow and external body force [195], but ignoring the complications arising from 
counterion screening in electrophoretic flows. We used the fluctuating LB model 
(Sec. 3) in conjunction with the point-force coupling scheme described in Sec. 4.5. 

We were surprised to find that the polymer migrates towards the channel center 
under the action of a body-force alone, while in combination with a pressure-driven 
flow the polymer can move either towards the channel wall or towards the channel 
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Fig. 1 1 . Center of mass distributions for countercurrent application of an external body force 
and pressure-driven flow in a channel of width H . The solid curve shows the level of migration 
under the pressure-driven flow only. The flow Peclet number in all cases is Pej = 12.5. The 
boundary is at y/H = and the center of the channel is at y/H = 0.5; H = 8R g . 

center. The external field and pressure gradient result in two different Peclet num- 
bers: Pe = URg/D and Pe f = ^R g 2 /D. Here U is the average polymer velocity 
with respect to the fluid, and 7 is the average shear rate. The interplay between force 
and flow can lead to a wide variety of steady-state distributions of the polymer center 
of mass across the channel [195]. For example, in a countercurrent application of the 
two fields, the polymer tends to orient in different quadrants depending on the rela- 
tive magnitude of the two driving forces. The polymer then drifts either towards the 
walls or towards the center depending on its mean orientation. The results in Fig. 1 1 
show migration towards the boundaries when the external force is small (Pe < 30), 
but increasing the force eventually reverses the orientation of the polymer and the 
polymer again migrates towards the center (Pe > 100). 

The simulations mimic recent experimental observations of the migration of 
DNA in combined electric and pressure-driven flow fields [ 1 96, 1 97] . The similarities 
between these results suggest that hydrodynamic interactions in polyelectrolyte so- 
lutions are only partially screened. In fact, within the Debye-Hiickel approximation, 
there is a residual dipolar flow field [198]. Although this flow is weak in comparison 
to the electrophoretic velocity, its dipolar orientation enables it to drive a transverse 
migration of the polymer. A recent kinetic theory calculation [199] supports these 
qualitative observations. 
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5.4 

Sedimentation 

The previous examples have focused on sub-micron sized particles, colloids and 
polymers, where Brownian motion is an essential component of the dynamics. For 
particles larger than a few microns, Brownian motion is negligible under normal 
laboratory conditions and a suspension of such particles can be simulated using the 
deterministic version of the LB equation (see Sec. 3). There is an interesting regime 
of particle sizes, from 1-100 micron depending on solvent, where Brownian motion 
is negligible, yet inertial effects are still unimportant. This corresponds to the region 
of low Reynolds number (Re = U d/rjkin) but high Peclet number (Pe = Ud/D). 
Here U is the characteristic particle velocity, d is the diameter, and D is the particle 
diffusion coefficient. Because of the large difference in time scale between diffu- 
sion of momentum and particle diffusion, it is quite feasible for Pe to be 6 — 10 
orders of magnitude larger than Re. We have carried out a number of simulations 
in this regime, with the aim of elucidating the role of suspension microstructure in 
controlling the amplitude of the velocity fluctuations as the suspension settles. 

In a sedimenting suspension, spatial and temporal variations in particle concen- 
tration drive large fluctuations in the particle velocities, of the same order of mag- 
nitude as the mean settling velocity. For particles larger than a few microns, this 
hydrodynamic diffusion dominates the thermal Brownian motion, and in the absence 
of inertia (Re <C 1), the particle velocities are determined entirely by the instanta- 
neous particle positions. If the particles are randomly distributed, then the velocity 
fluctuations will diverge with increasing container size [200], although the density 
fluctuations may eventually drain out of the system by convection [201]. However, 
experimental measurements indicate that the velocity fluctuations converge to a fi- 
nite value as the container dimensions are increased [202, 203], but the mechanism 
by which the velocity fluctuations saturate is not yet clear. Some time ago, Koch 
and Shaqfeh suggested that the distribution of pairs of particles could be modified 
by shearing forces induced by the motion of a third particle, and that these changes 
in microstructure could in turn lead to a screening of the long-range hydrodynamic 
interactions driving the velocity fluctuations [204]. However, detailed numerical sim- 
ulations found no evidence of the predicted microstructural changes [205]. Instead 
the velocity fluctuations in homogeneous suspensions (with periodic boundary con- 
ditions) were found to diverge with increasing cell size. More recently, it has been 
proposed that long wavelength density fluctuations can be suppressed by a convec- 
tion diffusion mechanism [206, 207], but a bulk screening mechanism cannot be 
reconciled with the results of computer simulations [144, 205]. Alternatively, it has 
been suggested that the vertical walls of the container may modify, although not 
eliminate, the divergence of the velocity fluctuations [208]. Most recently, it has 
been shown [209, 210] that a small vertical density gradient can damp out diverging 
velocity fluctuations. 

Lattice-Boltzmann simulations were used to test these theoretical ideas, com- 
paring the behavior of the velocity fluctuations in three different geometries [211]. 
We found striking differences in the level of velocity fluctuations, depending on the 
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Fig. 12. Relative velocity fluctuations (AU^ / (U z ) 2 as a function of container width. The 
vertical ((AU^) and horizontal ((AU%^) fluctuations are shown for three different boundary 
conditions: a cell bounded in all three directions by no-slip walls (Box), a cell bounded by 
vertical walls (Walls), and a cell that is periodic in all three directions (PBC) [144, 205]. The 
statistical errors are comparable to the size of the symbols. 



macroscopic boundary conditions. In a geometry similar to those used in laboratory 
experiments [202, 203], namely a rigid container bounded in all three directions, 
we found that the calculated velocity fluctuations saturate with increasing container 
dimensions, as observed experimentally, but contrary to earlier simulations with pe- 
riodic boundary conditions [144, 205]. The main result is illustrated in Fig. 12, and 
suggests that the velocity fluctuations in a bounded container are independent of 
container width for sufficiently large containers. On the other hand, in vertically ho- 
mogeneous suspensions velocity fluctuations are proportional to the container width, 
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regardless of the boundary conditions in the horizontal plane. The significance of 
this result is that it establishes that velocity fluctuations in a sedimenting suspen- 
sion depend on the macroscopic boundary conditions and that laboratory measure- 
ments [202, 203, 212] are not necessarily characteristic of a uniform suspension, as 
had been supposed. Instead, the simulations show that vertical variations in particle 
concentration are responsible for suppressing the velocity fluctuations, which oth- 
erwise diverge with increasing container size, in agreement with theory [200] and 
earlier simulations [144, 205]. The upper and lower boundaries apparently act as 
sinks for the fluctuation energy [201], while in homogeneous suspensions velocity 
fluctuations remain proportional to the system size [200]. 

5.5 

Inertial migration in pressure-driven flow 

At still larger particle sizes, typically in excess of 100 microns, the inertia of the 
fluid can no longer be ignored. For suspended particles in a gravitational field, the 
Reynolds number grows in proportion to the cube of the particle size. Inertia breaks 
the symmetry inherent in Stokes flow and leads to new phenomena, and in particular 
the possibility of lateral migration of particles. A particle in a shear flow experiences 
a transverse force at non-zero Re, with a direction that depends on the velocity of 
the particle with respect to the fluid velocity at its center. Thus if the particle is mov- 
ing slightly faster than the fluid it moves crosswise to the flow in the direction of 
lower fluid velocity and vice-versa [213]; if it is moving with the local stream veloc- 
ity then it does not migrate in the lateral direction at all. Now in Poiseuille flow, a 
spherical particle moves faster than the surrounding fluid because of the Faxen force, 
proportional to the curvature in the fluid velocity field. Thus particles tend to migrate 
towards the channel walls [214]. However, near the wall the particle is slowed down 
by the additional drag with the wall and so eventually migrates the other way. At 
small Reynolds numbers (Re < 100), these forces balance when the particle is at a 
radial position of roughly 0.6R, where R is the radius of the pipe. In a cylindrical 
pipe, a uniformly distributed suspension of particles rearranges to form a stable ring 
located at approximately 0.6R [215]. Theoretical calculations for small particles in 
plane Poiseuille flow give similar equilibrium positions to those observed experimen- 
tally [216, 217]. The profile of the lateral force across the channel shows only one 
equilibrium position, which shifts closer to the boundary wall as the Reynolds num- 
ber increases. Our interest in this problem was sparked by two recent experimental 
observations: first that particles tend to align near the walls to make linear chains of 
more or less equally-spaced particles [215, 218], and second that at high Reynolds 
numbers (Re ~ 1000) an additional inner ring of particles was observed when the 
ratio of particle diameter d to cylinder diameter D was of the order of 1:10 [219]. 
Large particles introduce an additional Reynolds number, Re p = Re ■ (d/D) 2 , which 
may not be small, as assumed theoretically [216, 217]. We used the LB method to 
investigate inertial migration of neutrally buoyant particles in the range of Reynolds 
numbers from 100 to 1000 [220]. Individual particles in a channel with a square 
cross section migrate to one of a small number of equilibrium positions in the cross- 
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sectional plane, located near an edge or at the center of a face; we could not identify 
any stable postions for single particles near the center of the channel. 




Fig. 13. Snapshots of particle configurations in a duct flow at different Reynolds numbers; 
the flow is into the plane of the paper: (a) Initial configuration (b) Re = 100 (c) Re — 500 (d) 
Re = 1000. The ratio H/d = 9.1, the number of particles N — 32 and the volume fraction 

<f> = 1%. 

To investigate multiparticle suspensions, random configurations of particles were 
prepared at a volume fraction = 1% and size ratio H/d = 9.1. The Reynolds num- 
ber in the simulations varied between 100 and 1000. An initially uniform distribution, 
shown in Fig. 13a, evolves into three different steady-state distributions depending 
on Re. At Re = 100 (Fig. 13b) particles are gathered around the eight equilibrium 
positions and strongly aligned in the direction of the flow, making linear chains of 
more or less uniformly spaced particles. Similar trains of particles were observed in 
laboratory experiments [218]. At Re = 500 (Fig. 13c) the particles are gathered in 
one of the four most stable positions, near each corner. By a Reynolds number of 
500, the trains are unstable and the spacing between the particles is no longer uni- 
form. Instead transient aggregates of closely-spaced particles are formed, again near 
the corners of the duct. However at still higher Reynolds number, Re — 1000, there 
is another change in particle configuration (Fig. 13d), and particles appear in the 
center of the duct. A central band was first observed in experiments in a cylindrical 
pipe [219], but its origin remains unclear. We observe that the central particles have a 
substantial diffusive motion in the velocity-gradient plane, whereas the particle trains 
exhibit little transverse diffusion. Since there are no single-particle equilibrium po- 
sitions at the duct center, the presence of particles in the inner region is clearly due 
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to multi-body interactions. Nevertheless this migration cannot be a shear-induced 
migration of the kind that occurs in low-Reynolds number flows [221]. 

Our simulations suggest that the inner band of particles is the result of the forma- 
tion of transient clusters of particles. We proposed [220] that at higher Reynolds 
numbers (Re > 500) the trains become unstable and clusters of closely-spaced 
particles arise, as can be seen in Fig. 13c. Simulations of tethered pairs of par- 
ticles have shown that additional equilibrium positions arise for pairs of particles 
at Reynolds numbers in excess of 750 [220]. Thus transient clusters are formed at 
higher Reynolds numbers, which drift towards the center of the channel making the 
additional ring observed in experiments [219] and simulations (Fig. 13d). Eventu- 
ally, the cluster disintegrates from hydrodynamic dispersion and the particles return 
to the walls. At steady state, there is a flux of pairs and triplets of particles moving 
towards the center, balanced by individual particles moving towards the walls. This 
also explains why the particles in the inner region are highly mobile, while those 
near the walls have a very small diffusivity. 
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